{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical Hypothesis Testing with Kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Mean Discrepancy: do X and Y have the same distribution?\n",
    "\n",
    "Let's suppose we are given two sets of samples $\\{X_1, \\ldots, X_N\\}$ and $\\{Y_1, \\ldots, Y_M\\}$ draw iid from distributions $\\mathbb{P}_X$ and $\\mathbb{P}_Y$ respectively. \n",
    "\n",
    "We'd like to know whether in fact $\\mathbb{P}_X$ and $\\mathbb{P}_Y$ are the same distribution. We can do this using the kernel mean embedding.\n",
    "\n",
    "The idea is this: embed the samples $X_n$ and $Y_m$ into an RKHS and look at their means $\\tilde{\\mu}_X$ and $\\tilde{\\mu}_Y$. As $N$ and $M$ become large, these means will converge to the mean embeddings of $\\mathbb{P}_X$ and $\\mathbb{P}_Y$ ($\\mu_X$ and $\\mu_Y$). Provided that we use a _characteristic_ kernel, the mean embeddings of $\\mu_X$ and $\\mu_Y$ will be equal if and only if the distributions are themselves equal. Of course, we only have access to finite samples, and so even if $\\mathbb{P}_X$ and $\\mathbb{P}_Y$ were equal, it would be very unlikely that the empirical embeddings $\\tilde{\\mu}_X$ and $\\tilde{\\mu}_Y$ would be equal, so we'll then have to do some clever sampling method to find the distribution of this difference under the hypothesis that $\\mathbb{P}_X$ and $\\mathbb{P}_Y$ are equal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a kernel $k$, the mean embeddings are defined as\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mu_X &= \\mathbb{E}_{X\\sim\\mathbb{P}_X} k(X, \\cdot) \\\\\n",
    "\\mu_Y &= \\mathbb{E}_{Y\\sim\\mathbb{P}_Y} k(Y, \\cdot) \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The empirical embeddings are\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\tilde{\\mu}_X &= \\frac{1}{N}\\sum_n k(X_n, \\cdot) \\\\\n",
    "\\tilde{\\mu}_Y &= \\frac{1}{M}\\sum_m k(Y_m, \\cdot) \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The RKHS distance between these two can actually be evaluated:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\|\\tilde{\\mu}_X - \\tilde{\\mu}_Y\\|^2_{\\mathcal{H}_k} &= \\langle \\tilde{\\mu}_X, \\tilde{\\mu}_X \\rangle - 2 \\langle \\tilde{\\mu}_X, \\tilde{\\mu}_Y \\rangle + \\langle \\tilde{\\mu}_Y, \\tilde{\\mu}_Y \\rangle \\\\\n",
    "&= \\frac{1}{N^2}\\sum_{n,n'} k(X_n, X_{n'}) - \\frac{2}{NM}\\sum_{n,m} k(X_n, Y_m)  + \\frac{1}{M^2}\\sum_{m,m'} k(Y_m, Y_{m'})\\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: this will be a biased estimator because of the fact that in the $XX$ and $YY$ terms there are $N$ and $M$ terms where the same datum appears in both the left and right argument."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So if we construct a big matrix $K$ with entries\n",
    "\n",
    "$$\n",
    "K = \n",
    "\\begin{pmatrix}\n",
    "K_{XX} & K_{YX} \\\\\n",
    "K_{XY} & K_{YY}\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "then our test statistic is just the averages of the diagonal blocks minus the averages of the off-diagonal blocks.\n",
    "\n",
    "This is known as the MMD (maximum mean discrepancy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def gaussian_kernel(Z, length):\n",
    "    # Z.shape = [N, K]\n",
    "    # broadcasts Z to be NxKx1 tensor, Z.T to be 1xKxN\n",
    "    Z = Z[:,:,None]\n",
    "    # then Z - Z.T is NxKxN, and we sum over the data dimension K (axis=1)\n",
    "    pre_exp = ((Z - Z.T)**2).sum(axis=1)\n",
    "    \n",
    "    \n",
    "    return np.exp(-(pre_exp / length))\n",
    "\n",
    "def MMD(X, Y, kernel, kernel_parameters):\n",
    "    N = len(X)\n",
    "    M = len(Y)\n",
    "    Z = np.concatenate([X, Y])\n",
    "    K = kernel(Z, kernel_parameters)\n",
    "\n",
    "    KXX = K[0:N, 0:N]\n",
    "    KXY = K[N:, 0:N]\n",
    "    KYY = K[N:, N:]\n",
    "    \n",
    "    return (KXX.sum()/(N**2))  - (2*KXY.sum()/(N*M)) + (KYY.sum()/(M**2))\n",
    "     \n",
    "\n",
    "def MMD_test(X, Y, kernel, kernel_parameters):\n",
    "    mmd = MMD(X, Y, kernel, kernel_parameters)\n",
    "    n_samples = 100\n",
    "    null_dist = MMD_null(X, Y, kernel, kernel_parameters, n_samples)\n",
    "    p_value = (null_dist[:,None] > mmd).sum() / float(n_samples)\n",
    "    \n",
    "    return mmd, p_value\n",
    "\n",
    "def MMD_null(X, Y, kernel, kernel_parameters, n_samples):\n",
    "    S = [False for _ in range(n_samples)]\n",
    "    Z = np.concatenate([X,Y])\n",
    "    for i in range(n_samples):\n",
    "        np.random.shuffle(Z)\n",
    "        S[i] = MMD(Z[0:len(X)], Z[len(X):], kernel, kernel_parameters)\n",
    "    S = np.array(S)    \n",
    "    return S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+cVOV1/z9nZu+ys5CwCKS4CwSTr8VGJFDR2rBNqkSx\nseBG45pE2+TbGtLaxsArhaBJYKEm8uPVajD1W/1imlhtAqhZl5h80SiNhZQquAtKIm1iirCLEYTF\nwsyyMzvP94+ZO3vn3ue597kzd37dOe/XyxfO7J07z9yZe57nOedzziEhBBiGYZjwEKn0ABiGYZhg\nYcPOMAwTMtiwMwzDhAw27AzDMCGDDTvDMEzIYMPOMAwTMtiwMwzDhAw27AzDMCGDDTvDMEzIaKjE\nm06aNEnMmDGjEm/NMAxTs+zbt++EEGKy13EVMewzZszA3r17K/HWDMMwNQsRHdY5jl0xDMMwIYMN\nO8MwTMhgw84wDBMy2LAzDMOEDDbsDMMwIYMNO8MwTMioiNyRYZjapru3Hxt3HMLAYAKtLTEsXzgT\nHXPbKj0sJkvRhp2ImgC8AGBM9nyPCyFWF3tehmGqk+7eftz55CtIJEcAAP2DCdz55CsAwMa9SgjC\nFXMOwFVCiA8CmAPgWiK6IoDzMgxThWzccShn1E0SyRFs3HGoQiNi7BS9YheZbthnsg+N7H/cIZth\nQsrAYMLX80z5CSR4SkRRIuoD8BaAZ4UQ/xHEeRmGqT5aW2K+nmfKTyCGXQgxIoSYA2AqgMuJaJb9\nGCJaQkR7iWjv8ePHg3hbhmEqwPKFMxEzonnPxYwoli+cWaERMXYClTsKIQYB7ARwreRvDwkh5gkh\n5k2e7FmcjGGYKqVjbhvuueEStLXEQADaWmK454ZLOHBaRQShipkMICmEGCSiGICrAawvemQMU2JY\nslc4HXPb+FpVMUHo2M8H8F0iiiKzA9gqhPhhAOdlmJLBkj0mzAShijkAYG4AY2GYsuEm2WPDztQ6\nXFKAqUtYsseEGTbsTF3Ckj0mzLBhZ+oSluwxYYaLgDF1ielHZ1UME0bYsDN1SxCSPZZMMtUIG3aG\nKRCWTDLVCht2himQMEsmu3v7sWb7QZyKJwEALTEDXYsvrvnPVS+wYWeYAgmrZLK7tx/LH9+P5Mho\nkdbBRBLLt+0HwLuRWoBVMQxTIGGVTG7ccSjPqJsk02K05vqBrcC9s4Culsy/B7aWeZSMG2zYGaZA\nwiqZdNtxDAwmMkZ8+x3A6SMARObf7Xewca8i2LAzTIGEtcqh246jtSUGPLcWSNqMfzKReZ6pCtjH\nzjBFYJdMdvf2Y/6652ta/rh84UyHjx0AjAhldiNPHZW/8LTieabssGFnGBuFatMrJX8MWktvvlap\nivnXqVk3jI3xUwt+TyZY2LAzjIVijLNXk+dSJDKVajJxTd5asCrjU7e6Y4xY5nmmKmAfO8NY8DLO\nbqiCjqax7R9MQFged/f2V3S8BTO7E1i0CRg/DQBl/l20KfM8UxXwip1hLBSjTW9tiaFfclyUqGSJ\nTBXT0s/uZENexfCKnWEsFKNNV8kfR4RTEw4EY3zDqqVnioMNO1PTmCqUC1Y+jfnrni/avVGMNt0q\nf7w+sgt7mr6In0c/iT1NX8TiyC7H8RGiio6XCS/simFqllIEDost59sxtw0d0d3A9n/KBRen4DjW\nG5uBJNCTbs8dOyIE7nzyFew9fBI7Xzte8PsVM14mnJBQbBNLybx588TevXvL/r5MuJi/7nmpT7ut\nJYbdK6+qwIiy3DtLKgfsF5Mw/9wmx/MEwHoXxoxoKBKdmOAhon1CiHlex7ErhqlZqrYIlyJR53y8\nLX3evrQquaqFCT1FG3YimkZEO4no50R0kIi+GMTAGMaLqg0cKhJ13qJJ2qeo+OTE1DRBrNhTAL4k\nhPgAgCsA/BURfSCA8zKMK1UbOFywKpOwY8WI4cjvLneMlxSnqPjkxNQ0RRt2IcQxIcTL2f//HwC/\nAMDOQabkVG0RLkUCz2WLP+8Y7y1XTHcY+080/gzP0u1cEpcpmECDp0Q0A8ALAGYJId5RHcfBU4YZ\nxVrr5TPjXsRXxT+iYWRo9AAjVtbMzqJrzxzYmqn0ePpoxi21YBUnMwWEbvA0MMNOROMA/BTA14UQ\nT0r+vgTAEgCYPn36pYcPHw7kfRkGCFFTaYWiBuOnActeLfnb2yWkgE+Vjlmr3V5HhksOBEJZVTFE\nZAB4AsBjMqMOAEKIh4QQ84QQ8yZPnhzE2zIMgFFjVIpaLGVHVfq2TCVxi649w7Xaq4IgVDEE4GEA\nvxBC/H3xQ2IYf5SiEFbQGa3aqErflqkkbtES0gpPTEyGIFbs8wH8CYCriKgv+9/HAjgvw2gRtJ69\nojsAhaLGrSRukJNQ0RLSCk9MTIaiSwoIIXZBrdpiQkY1+rJVVRULlQy67QCC/Kzya5n1Q2sGH4Mu\nq7B84Uypj11bQsq12qsCrhVTDqpMJVBrHYK8KMQYuV2DQnYAfq+p+7XUL4kb9CRUdO2Z2f4mJqY0\nsGEvNXaVgNnRHajIj71UHYIqadj9GiOva+B3B1DINQ3qWpairIJr9yQduFZ7xWHDXmrcVAIV+PEX\nY1CqtjYL/Bkj1TXo6jmIjTsOoX8wIS3MpdoBFHJN/VxLt92AahIaHzNqvqk2UzhcBKzUVJlKoNgO\nQX6er1ZUn3UwkcwZSYHRwJFXRmsh11T3WnoFcmVlFYwI4exwKhzyT6Yg2LCXmipTCZSiQ1DFa7P4\npKXZ0DpOYLQEcMfcNqX6pJBrqnstvaScsrIK45oakBwRytcw4YddMaWmylQCyxfOxPJt+5FMj974\nRoS0OwQBtd3Uobu3H2eGUtrHm6tuNz96IcFb3Wupsxuwu6EuWPm052uCoBoVUkwGNuylphpVAnZx\nqg+xatGBtQrT1XMwb1IzIQJk1TXMVbfbytls6uHXyOlcy0KknC3NBk7Fk75e45dqVUgxGdiwl4Mq\nUgls3HHIsU1PjoiKK1vKQXdvPwYTToMHZIx6zIgqV91eK+dSTXh+dwOqHYkR1duV6VKtCikmAxv2\nOqMalS1eW/qgtvxuPua27Hn9qk9KHTh2c9nIrsvGHYekO5KxjQ2BGtxq/B0xo7BhrzMK0WgXY1R1\njLZ9S79z27fwkR89gQnJtxCPTcGuszeif/hDub8XuuV3MzrmuFTn1F05l8LvLBuXyhViX0WbnFbs\nVAqlUhMdowerYuoMP8qWYmum6LzevqVfHNmFe4zNmJD8DQCB5sQxrKWHsDiyK3dMoQoPldGZ0Gxo\n+cO9mnqUs8aMyhUSJXnAJGiDGxaFVFgJtNGGLtxoo7KoVpX25+PDKWkQzpQAejF/3fPSVZ319Res\nfDovEWhX4x2YGjnheM3R9CS0D2/KPSYAv153nfZnM/9mX9UaUcLYxgacTiSLXmHrfN6gsF83K9ZY\nweLILnzZ2IpWehsUcOCeVTHlR7ceO7ti6hDdrb0KXT+q6jjrue1b+lZyGvXM82/nP5asQL2UGnZ/\ndUuzgTNDqVxA1Tx+7+GT2Pnacd8Gq5x+Z5UrxBormPfOs1jX+DBiOJf5Y8DlLKpFIcUTjBN2xSio\nWD3uCiHb2qtwGNUDWzOdf2w9OlXbfwLyMietzoMBMUn6mgExMff/qi2/Tl32jrlt2L3yKvx63XVo\nbmxwBBoTyRE8tueNgtwp42PyxCfdhCg/uLlCzM/4zcnbR426SciaXoSqyUqAsGGXUI8/Ft1VpcOo\nmkXOTh8BIEZXhQe2Ooy2iQCwZvtBzF/3PJZt6UOTMfoz3JDqRFw05h2fijZhc+Otng2r/a6YVc/b\nXRw6Pv3u3n6cHZYnPp0ZSgX+29Fq5F1l5SxKQSmarIQBdsVIqEeNrmpr3xIzMHZMg3qb61LkrGPZ\nq1i6pU/6fqfiyZz/PpFMw4gQxjU1YHu8HecZjVhhbEFz4k1g/FQ0LFiFrtmd6CrwM8jqr2zccUjp\no5bhNfGt2X7QkR9gkkyLXIGxUqtlAIyWiVZ9whKXsyina4Rll3LYsEsI649FdsMBcK1o2LX4Yveb\n0mNV2KYwtnaSaYHmxgb0rroGwHUA1mh9jvExA0TAYDyJlmYDRoTy3Cv2HYYsgKqDm6rkpZ4HsT21\nAa1jTmBATMKGVCd60u15xwwmkg5fPlCCLE1ZM2krJS5nUe6MVJZdymFXjISwVDG0InMvLd+2H8sf\n3y+taNgSM9BkRLB0Sx/ef+ePMEMVa/AocibzBavQmTjtn2MwkVn5C2R2AaDM2FXuCT+xBBNXGd+B\nrZj18tcwNXICEQKmRk5gnbE5T54po2TuAtkOymT8NGDRppJmQZfbNcKySzm8YpdQdHuwKkR2w8ky\nFAUyhvFcKp07fiQriZWuvjyKnJnHdfUczK1Yveqy+P0ceZ9pRGDsmAb0rb5G+ne/u64oEW681EX9\n8dxaR4CymYaxomEreobb5a8pcCyAhptD6T8nYNmrvt/PL+Xe7YahMF0pYMMuIYw/Fj83lqqeCiCJ\nNbgUOTONkOnmMZEZdbckKev3oOPW8aqDrnMOkxEh8MS+fsx773ny719hSO3yTNVYvLC7nc4Op3K+\nfOlEO35qNpBtw7azKpUfvBKukWqRXVYTbNgVhO3H4tegueEwnPYiZwe2Ir7+IiyOv4l5YiI2RJw+\nZyCzGk4LIU0kkvn9ZXEAGV510P362M3uStLfg8KQWuWZMnR2gHZ/tWzCdUy0GmWiS+kHD+NutxYJ\nxMdORN8moreIqPR7PaYgVJ12jGi+IDFmRDHBQ3ftVlem6+7ViD/xV2hOHEOEhKvPOS0Efr3uulwj\nC/Mcpg8dcBpxL6MuC5Za8xEA5MkEJ2QDrl4MJpJyyeKCVRnDaSGBRmxIqf3YE5qNjO8/uluq/zfR\njQfkTbSzOzN+9PHTAJDUr15KP7iWDJMpOUGt2L8D4FsAHgnofEwJaDIiuRu6JWaga/HFAJwup72H\nT+KxPW9IjahXXZln6VE0R4bz/qbyOcsmCF1jJlu5T2g2sHrRxY5Jwr4yveeGS/JS/Lt7+/NiACqk\nq3aJK2rt6Y+jJ32F8jzNjQ0Zo+7R5Fx3h+W4jh5lokvtBw/bbrcWCcSwCyFeIKIZQZyLCR6ZxO9c\nKg3AeRN29/bjiX39UqPe5uKLNQ1y6xi9kgCqCULHuHiV2LWPyUoiOYI1250G2rwebpirdqlxtxjS\n7yk6GJkMDCY8m5x39/ZruZ0KcXOU0w/O6f6VoWw+diJaAmAJAEyfPr1cb8vAX8KVasXsVcjKNBQD\nYhKmSuq9DIiJOUPlNkF4xQKsafNeBkI1SZyK5xtoPxLIpVv6sHHHoaIMVGtLzFP/75ZANaHZwGC8\n8KJl5fKDc5elylE2wy6EeAjAQ0CmumO53pfxt/UuZJtuXV1uSHVinbEZzTTqjomLRmxuvBX3fnyO\n5w0tMzo6E4IMt0nCOqn5dUGoDJS5OpWxOLILKxq2opVOYIjOB2ITgMRJ54FZ9YrbmDJJXIVTLtVX\nPWZwVwusiqkD/Gy9laUFXAKq1tVlT7odSCJrxN7GUPMUNP/RWswZmY+NOw5h2ZY+V0MSpNFZvnCm\nsqSB1XBaP7PVAKuySAGngXLLaF0c2ZU32TUnjgERA4g2AiP5E+CGszdiTm+/a/VGGX5dHqXyg1vH\noVq91XoGdy3Amad1gJ/svOULZzqUMoB7ISv7jdqTbkf78Ca879xjaP7ya+geme+rqJq1AqPp/imk\n0mbH3Da0KCouWic18/qYBlg3i9T6ud3cOXc1bsvbwQAA0kmgcRzisfORFoSj6UnYNvJh3Db8KBZ3\nX4xnkp/Fy2OW4PUxn8auxjuwOLKrZA1RgsI+DhW1nMFdKwQld/wegH8HMJOIjhLRnwdxXiYY/EjQ\nOua2YWyjcyOXTAulm0F1o5qry2LkdTKjtWxLH77a/YrnawGga/HFnpOaeX1kBthU9MiIEOWMp2wV\nujiyC7sa78AUHJcPLnEK889lJsANqU7cFH0hO6kIjB05jfPoTG6CWd/4MB657HCuIYp1ouvqOVjS\nNH7dEtY6sQrWtJeHoFQxnwriPEzp8LP1VvXHVG2hvYJxbg037CoTu0vh7LmUw1gIAI/teUOdDWo7\nl9kybkQIpZ++Y24b8JSeosdkRIicr93uOrG7X2TEY1Nw6lTmWq9o2Op6bAzncNmv7kd37x8H0hBF\n13XjJwDq9n4EsCqmjLArhnHgtwia147AbettdRnIVuemttxc/ZquiUWRXa4rUnui04gQeYoa2fFv\nwrvJhx1zZWx3d3kZ6lS0CRuSN+ceqzpH5XH6qC8FjwCkK2w/rhs/uy23nZs9EY0pLWzYGQeFVMzr\nmNuG5QtnorUlhoHBBDbuOJTXJUlV4dFqJEwjYjfiaxq+LfV9z3vnWeV4VAbpS1v3O1wKpqH7xvBN\njiYfCeGeRQpkVqrm5GY2k1YZaiGAlIggkhrK+NOz/ntV5ygraSJc6vKZZeg0EAdGvwe720W1I5Ct\nzrnSYvXAzawZKboNr6+8aDJ2vnZcWc/dXLl39/YrFSoA8N/rrsMFK5/GIokLIy0AWdb/m5iMPdf/\nVDpOt2bP9vGZdWkAqyrmbQyIidiQ6sR20Q4hkHPn2JE151Y15bZ/lrhoxMrkbQDg6bqxHi9T6gSB\ntRE2IM/wBdR5DZyQVFp0m1mzYWe0KaRJhdUAqFaABODem+dg445D2BL/nNQgyhAgfGDk+w7fvt1Y\nu6Ey1iq8Jq8vbd2PESGkPnbVBHU0PQntw5vyJpVBjMV4nIFEoJQ7PmhU18LtMzPlhQ07Uxxme7XT\nRzPJNADSiVMYSE9UartlEIBfr7sOQGZiWLalT7oCNI3K62M+LTV+MgbxLpxJj5FqzpuNCOJJdZkA\nXb266jMJAJ8d92KuhV88NgWrzt6Ix4c/JHmPzOq/lU5IP1sahPcNPeZ4XnUt0oLwvnPO461jk+H2\nme0rdfvxb9Fk3DN8E/a++2pehVcQXcPOPnbGib1BdeIkkDiJCNyrNcqwBtQ65rYpjY65UlT5mh09\nQSIGmkVCqTn3Mup+9Op2RPYcK5IPZJKNINCcOIa19FDeOax6/vbhTcrPdkwRnD2mEcxtiRl5Qetb\nrpgujWfIPvP6xodxfWRXLthtTX6yHz8Fx/HNsf+E3R87wUa9Bqgbw66rxa1H7Ncm/uNV6vZqcNd2\nW5EFzszgoooNqU5HADMuGvHPIx/Fm5iMXCnaMe9CI6UKGpdMsaL7WiBj9P7e+Eff55B9tgQasT7p\nDM5GiTBw6QpHSeC4JZhrRAhdiy/OS+a6u+OSPIWSiewzx3AO35y8PadWsQY/paoes0gZU/XURUkB\nLkakRnZtmsa8CXi4Q1Tabq+6Ll7+bHtJAjOA2ZNux+pU9rx/OBMdT13sa1z5x+jr1e3ui+fSc3BT\n9AU0kHxH4Pb+bp/NTloIXLb488CMCcBzayFOH8WAmJibBHY13oFWOoHTP/otdD39KXz3zOV5wUrz\nut/yf/8du391Ui2ntBQjM1+zZvtBtI54H18qOABbPHVh2LkYUQbZDSO7NgNiorRCo/2YNosqRvcm\nbNPo5NSTblf2CzUn5WvGTcm6QZzj8sKtAqUVewB0Kp3An9BPXGMA1nNcH9mF5RIj7tULFbC4sLIl\ngduzgWf7mCYkf4MV4gGcjAyjZ7Adu37wAK555omc3//C5M3YjcuVn1nWjPzMuRQGovrHBwkvwoKh\nLgx7uRvsViOqG0amcJFVaMzDiGHqonuwe7a6jK+KQlrT2UkkR7AheTO6jAfzXEZxDc05IP98aQE8\nl54DIFMW91Q8KXVHuBn1VLQJmxtuBQ0Dnxn3Ir4qvo2GkSEAmUlhnbEZSCJvhW5vHA6MurBkBbVc\n3UgpYC1tRnNitNiYafQ3pDqx3tiMmPW1tpZ5QGa1nhwR2CAkvwHJ8XaKXW1XahEWtl1CXfjY/WZS\nhhHVDSPzefek27HBuH20vVrsvMx/ilZrfpBlqd56xfTcYy8fvMl3z1ye1wLuTUzW1nf3pNuxbeTD\neQHZCAE3RV/IBT/bWmJ62aBZ0hRBw/X3Y851S9DaEsNtw4/mjLqJ3QcfM6LoWnyxNGt37+GTWLal\nz1FQy82N5Gb0e9Lt+HLytvw4heR7PBVP5q7RyuRtOJqelCtS9tIla1y/9yCKkVViEVYtRdSCpC5W\n7LXeYDeI1YTqxjBT7e3XZs51S4C5a4oatwq3ujUXeHQfMmltiaF7ZA42ntuEgSH3aoLz338eXn7j\ndN5nXBDpc6y+c0Yw3o7rZp+PgZfl7ghZktGGhtuRev138MS+V7Q7SVm14PZ6OarWhCqXyikxFm0e\nsYOedDu2D7Xn5Kde2N1GsZeiuGeapINUliBW2+Xs7mQSRldtXazYa7nBblCrCbc6HuW8Nl7qJN0b\neDA+jOXb9nuWiJ3QbOCxz/2+Q87nFUDd+dpx/GzG7UqFjnUluzJ5G75z5nI8tueNnIFQSRtNH3xb\nS8xxjc1rs1Sh9QfkyppzIop30RBUmx2r39/r+sYMtUkw2wqqvr8gVtuVKEsQRldtXazYgdptsFvs\nasJc7atS/nXbzAVBd28/lm/bj2TWB9I/mMDybfsBjK5Yly+cqUxisnJ22NtHHzOiWL0oo57piO5G\nx5i1SDcdxUB6Ik6JcZhIZxyvMY1g/2ACXx68CC81LsEXxfedCh3J+1nHrOoktSHVKTVUulm9MmVN\nMw3hPMlnsb4nkFEsuRnI7t5+pBwJA/mciidz7hp7YDOI1Xa5ujtZqcQuodSEasUeRq16MasJe4VD\n6y1biV1LV8/BnFE3SaYFunoO5h67JTH5IUo0+vksCVdmktW7aAjnRP7K0BpABTLXa9vwh9A+vAlL\nk38JALjPeCDX+MKNPB91NgZwZ/I27Hv31dLr7qdqY0+6HR9J3p9LfGqB3KgLgVzcgQDccsV01+97\n445DSI4IRxE2t89qLeJ25UWTpceonjex37cAsHvlVbj35sx3sWxLX0nv5zAWLwvNij2sMqliVhN9\nTz+EZ+lRtI5xppCfPZfCsgAaM7thjw0MKuq8557PljF4fcwR32n+dtJC5Oq3XPHUXZiC/GvYSCmc\nEU0wMJJb3ZgB1H3p3857X5nsUaZwsdOTbsePU3+AjTd9EB1z2/BNl/H62fYT8vMBVH7339BkbE+3\ne/aKte7qCvms5th3viZvKKJ63nxv+327bEsflm7py9thlvJ+rsQuodTUlGF3CyKGMQACFBH4PbA1\nk/Iekd+gpjEt1Q0ju2FNZDVLcOBsZlWdTGRS3jWNp4rWllhuDAcjx6UJV+NoyPGcVTpojjGNiCMh\nKRdo9dCkNzZEpFUx7YbDrfG2HTMBzDxeKk81Ypiy6Bv49Wz3QKn9e3JT1mwfbpfupsxFRiG7S9l9\nK2z/mpTyfq5VV62KmjHsXivyMAZAgCJWE8+tVUvfbMaoFDeMyrUgWxGuNzYDPx7nKGPgZjytk8Mp\nMQ5EQAvO5CaKqxb+NdZsz7SMG2iUr2gF5Am2rdlJxRxjBP6zTE3ODo9o7Sb96PvNxDBTOWP1u7dF\n3gaNn4qX3v8FLP3RJAz8y9PS34x1lW7/7KrPKjPqRoRyi4xCdpd+789av5/LRc0Ydq8VuepH1dJs\nYP6652t6i1XQakKR+q0yRn5uGB35pep80polNAyROKkwsvI0f6vhtQZBp9IJrG/cjFePzMCpeCZL\nUraiVRl1AEgj4lkXHdDLcgX0dpPmv109B5UuKyCzW7vyosl4Yl9/nqHdnm7Hu3/307i74xLLRJL5\nDuwTiVugVjcr12RcU4Pr5OS1u/SzUzGPZ7ypmeCp14pcFgAxooQzQ6lQJR5oo0j9Vt2gXjeMGeCa\nsfLpvMQZ1TVVJoX5SPoxxzuh2UBLzMg9t7rhEY/+oMOY9vLG3GN7sk1KRJRG/ZyIIqpYoVvRzXIl\nqPuS2n/THXPb0Lf6Gtx6xfRcohZRpgSxVYq687XjUvfFD/cfy8kl3drZuQVqVUXYVJ91MD46CVll\nxUAmgG3tzCTDrbuWnVoPaJaTQAw7EV1LRIeI6JdEtDKIc9rxyh6VadXHNjY4VBhBdm+vahasklYG\n/IfIp2HYMnO8bhg3dQ0gv6aqG1al7z6ZHucYL4wYpn7iHvSuugZ9q68BIbNaV0n7rLxH5E8g1hK6\nKteKEACBlHpw63HbRj6s5ft32xnIftPdvf14Yl9/LjgqBPAx2oWDE76E3UM3oONfFypbAg4mkq6r\nX/NvbrszcxI8Hn0P0sjo9DcYt+OFMVdKj48Q5anQrFUizc/gtqCyTwb2a2U+rqXck2qgaFcMEUUB\n/AOAqwEcBfASEfUIIX5e7Lmt6Gzz7C4LVRZjXfjpzNRvs1nG+KloXrAKvzcyHz+2bPcnNBtYvehi\nTxmcl+9Xtvo0X2s1Nip999+OfAb3LZqTN14sWJWXwt7aEsOK+FZPwwuM7kxkgVqVu2EEEUcpYBlE\nmcxVmZZdhtQ3HSXpZNrVczDvWi+O7Mqr/4LTR7CucTPEsP+gMiEzcXi5P3rS7dhnXI3dX7sKUwF0\nAZijcN/YjTfg3m922ZY+jI8ZIMqs9k1Xntllq5Q1W8JWD8aNIHzslwP4pRDidQAgou8DuB5AoIa9\nkCBiGBMPfJGtDGgi860OuTSkMNGZCGXX1Jxorf1H1aVr5+O+2dc5xrvREh+58qLJaO31DliargOV\ndO/xkQ/jE/SCY3KJwduvnvu8GoFTO9ZJ5hgmYeDICmDu53N/7+7td/jX5XXU9RQ5dgSQk7d6BWrd\nJuqBwQQikjZ65s7NrXwFgLzPaPf/l0qdElY5tIogXDFtAI5YHh/NPhc4HXPb8poKeH0hYUw8KAa3\nIJ4bOhOh7Jqafnn7itXeWagnq7W2v9ZeSuGJff043fge6fuPgPJS/HvS7Urp3lWRPkeBq5XJ29Cv\ncBPJIAitRCUTe0eiNjqBWS9/LaPdzyL7HlQxiTY6of3eVkydeJMRyYtbON5XMVGb919aUVffnIT9\noOseLSYBUfXb73v6IeDeWUBXS+bfA3rNVqqdsgVPiWgJEe0lor3Hj6sTFoKklmvElIJCJaFeAa4J\nzYa07okq1yWaAAAgAElEQVTVL++GVTJnoroRvyk+JY0d/G3DF3FpdGtuogDcpXuyyUUWOFT1BaFs\ne7lvGg9gX+MSLI7scu1Noupg9OaTd+UMlexaqWISRPDVzs+KQKY0wLlUGrdKWunpLH7cYl5+AqIm\nXr/BYmsmyc5vtjfMtYA8fSSTSxEC4x6EYe8HMM3yeGr2uTyEEA8JIeYJIeZNnuyeYhwkflf5YabQ\n8sUdc9tw46Wj182ecn5t+t8cN5jdV+xGMi2w1JY2rrrRzXK98dj5jiJcp+LJPOPqVYjLzr81XYmD\nl94NjJ+WCxw+MvJRvJ0eByHkRp4ImBg5g3XGZnyi8WeOnYeJapJ5jziRM1TWsZvXuJVOOPu9ZvHT\nzk9GIjmCna8dL2jx47Ybti+odEoxe/0GC91tup0/zO3/gvCxvwTgQiK6ABmD/kkAnw7gvEzAqALQ\nV1402VPrb6aFy/zWXxP/iFU/SAO4PaeTdtNiq7D6PV3jI7Ovw9U/moT+c86/6xbiknEqnsSSvguw\netEOLN3Sl3t+Nf4MAPD6mE8rV+XNNIwviu9j/sLbHdfYiBCOYRLa4K4PNxU0i2zXGMgqdiRvXoiv\nP+/9BxMF+bW9Yl7Wc3oVONPZIRSbgCj77SuvXRna/5Waog27ECJFRH8NYAeAKIBvCyEOeryMqQCy\nm9FMdvEKKs1751lsadyKNjrhMDDNNIx1eAB3P92AjrlripKTmqswLxWUzg3tp8eoyal4Enc++QqI\nnCt0ZXu5LK30ttLgDRxZgfNe/hpiOJc7XjbJCAB3NW5Dsy2Qq1OStxCKERLoTgj2ayJTxchcedZr\n2JLtalXo+GXfyxDJ2yuWuv1fOSDh0Vy4FMybN0/s3bu37O9bL/iRdal8u20tsZwEDQe2IvHkX+cZ\nJRlx0YjmG/8BF/zL2KIqNBKAX6+7Dt29/Viz/WDeDW3KM2Xp8Lp0RHfjb6Jb8iSQOtJB+27FzgAm\nobXrV+oTZIuc4fRRvIlJ+MbwTY73JQCvN90C0riCw6IBf5Ncoi17tJdtNis+znvveVUlA5St8I0I\nAQQkR0Y/QcyIFhczM6t+WktZGLGiOoSVGiLaJ4SY53VczZQUYPTwK+tSrXz7BxM598y/N92FKR5G\nHcis3PHcWrS2bCrY6AL5q7AzQ/m68lPxJJY/vh83XzZN2WXIjY837MbGMQ9Le5ECcGjerUazJ92O\nSIrw9aZH0Zw6nbeKTgrCxMZURl0h0eADyJOf7untx3aLu8dEAPgNJmEKvAUG/yOatI16zIjid6eP\nx89+dTKvyNaWl45gy4tH8mrkFyIDDFIjLvOnJ9MCLTEDY8c0BDcBSXI9pN9bDcKGPWT4rXKp8mVb\nU+HfI+TVEWWI00dxVpLkY0QJYxsbMJhISt0cJlZ3i6x+O5BZtT265w3dIQEY1ZC30QmQzdXbTMPo\nMh5BE4adxcls1SW7R+ajIUlYF30QDWJ0J9FAApQczDww1RXAqJGwrNYxfio6FqzCUoyVjvWe4Zvw\nzbH/5CiKZmcCndX67GbZ3o07DjkmQusK2MRvUTi/iwmvSUC12DidSKJv9TVaY9LGlusRFtiwhwy/\nQSaZL9u+ZffyLee9j5iIwaF8X6hbdmt3bz/6nn4Itw0/itbI2xiKTUFzdC26e+d7BmB1V+teLhQA\nmIAzDj92jIbxZWMres7lr4qX4vt5Rh2QzHvJBOI/XoXm2Z3OLX/W8H923OfxnTOXO8ay991XAx+b\nCzy3FuL0UYwIcpQNBvT866ZWXacrVd65fey4/CwmdCaBuk8sDICaKQJWDwTRAcqvpFGm9RfIlzTG\nMKRcYVtJYAzWJ52rn+bGBuXqryO6G130YCZxByITzNp+RyZxJCCksjZNzodTOaFbyKwp/mbmO3xu\nrXP1nUxghbFFrSGf3QksexXUNYjeS9chgTF5x9kDr0aUHDWAjAjh7PBoETw/+DGifhYTOrJFP4mF\nYeyaFgRs2KuEoJpWF5Jta9f6f3bci3lZkhMj6qJbQgBpQegXk/Dl4T+X+nxdV38Ko/eF4c3a7dm8\n8DLEcdGIk2Kc9G/H4FwVq/TxzuMmZgyWQj7XnHhTS0N+2eLPI3bDt/AmJjuya4GMYsYseGdqxtta\nYhjX1CB1tXjhNzvbz2JCZxLQTSwM6p4JI+yKqRKC6gAVRJuvFcYWNKf05Hb9YhLahzcBMBNRnIbE\ndfWnMHrnRc6Asr08vbopEYAmI6rUSatcSUJkxm+ufGWad9kORFrf3aYzN1fUA4MJ4LemZrMbbYyf\nqq8hn92J31eojYQYrb8yIkTOMC+TBGd18Ks08VOHXdfNonNdwto1LQjYsFcJuttZHfWB103hdY7m\nxJvS19lL0NrdAaZR8dXGb7zc6NnnEa9WdFMnNOG/3pIHE1WJStZVLwBtzbtMH/9ceg4WRPqkr+06\neyO+Gv3HnBIHQEZWt2BV3nmDap9nGjfV8WZ2rLQxTcxZHsILP4uJgls9Sghr17QgYB17laCjJ5fp\ne/1qed3OAWRuzi3xz2FqxLnCPWeMx4lhA+dDbvis6gtrAtTO1447b/icSuQInOFaOWlBeN+5x7Q+\np53Ryop6iUpB8/GG3fj6u57MTJoSWZ3Od+uVwWmFANx78xzX73r5tv0O1ZERJWz8xAdLuuINShqp\nlYMRMnR17GzYS4VN3ualj9W5sYP4IavOMaHZwFAyjURyRKoiiYtGbDBuR2rWJ/DD/cccihXZBKP6\nTGsuOIhFb6zPS3hKZ7cDQ7Hz0YxzQOKkY4xH06Nun6AwIpnys6p6LEHSEjOUcr05a56RqoDs363d\nKMaHU9KMTPN1bkb04lX/D2eHnZOEX8NYqTrnQSx0ag1OUKokCnkbAKVx19nOBrH1VB1rNQ7KVPxz\nlyO2rz9vdT8wmEBLswEhMpI6sxxAx9w2pQ/0Q//9AGKR/ISnCGUM99Vn7sMjlx3GZa+szguq6rai\n8wMRkAbKYtQBKOWbbrV1ZHXRdSZP07Whcst19/ZLjTowmpymY6Bl8sWlW/rQ1XMQXYvdG7gUSxDx\npLDCK/ZScO8sRbBsGrDs1YJPW8oVu4ms45DdZeHlHrImI8l4fcynEZEEY01XS1tLDLs/dgJ4bi3S\ng0cDc53YG13cJz6JbcMfKuqcfvnvddc5nnP7TnS+20JWzF6/A8DporNO4qcTSdcdg/X1bGiDg1fs\nlURVHa7IqnFBBJ5U5xjTEMGHz+2Udhyyq1Gsq0hp+veIcE0uUqlUzISbgcFELiPwDzwM0IRmAx84\n/13Y/Sun68bKxxt24xsNm3NdktpwAuvxLawf8y3PejE6k50OYxvlNcrddlw6320h1Rl1dnmJ5Ai6\neg7iXCqd+46tRtxrYghEoXJgK+I/XoWmxJsYSE/E5sZbMee6JTxZeMCGvRQolB6yqnF+Vlsd0d24\nZlxxP3LV9hUALuv+S2nHIbsaxSpN86tAIHiX07We/8qLJjtqwpjFq+7uuCT33AxFf1uTVbHHEUvm\nfzZz1+Amp1S111NJL90wos60ke7efmmbOWC0gUkpfNi6CptCyi9bKUqhcmArUk99Ac1mXZ/ICaxI\nPoBVP0jBLBHNyOEEpVKwYJWjy49K3qadYJH12zcnjiECgamRE+iiB9ER3e17eLLmIx1z25T1qa3P\n23cIfjIU21piuOWK6Xg2+hFpa7qedHve+bt7+/HEvn6HXkYAeGJff951UjW4MBk//BvXv6uaVqja\n6xXS4OK0zUia37/MqMeMKFYvurhkSTiFdDkqhKLKADy3Nl8iisy1X4rvF1Uauh7gFXsp0Kwa5yvB\nQpGhiefWBlbEiBQ7jbdoEgiQrhZ1GiMD+b7iTJnYRvzBYPtobe7hZE4uad1VqM5rv07LF87E8sf3\nKzMtderdyCY2t/Z6fhlv6zGq+nxRopxveu7aZ4pKwpGt9q3vHc3uFlpiBs4OpxxlcZuMiNKHbtKS\n/Q7txxXdX1jhumylt7V3ApVS7FQaNuylQqNqnC+VS9B+e5kcc8EqaX3qKYu+gV/Pdgb9gHzXjmpr\nb+9p6uYTNmt/DGjUN7GnoQNw1G83kbl/HOeTFNXyigeY6Cjx3xlKYs6aZ3KBR9X1SguRc8GojKqO\nYZMpVpZv259X13xECBCAP/7g+dK67IB7AbGYEc2pX1STiFd3LiWKhcaAmKi1E/BbdTJMsCumgvgq\n2KXo6hKPTfH/xqYc097EF8g0GRg/DQBl/tVoOmC6dlTukHFN6iJgVuxuBy9kaei9q66RKk960u05\n90+mvo3tACOGzY23Ol63IdWJBPIbXNullzEjiluumO453nQ29d90qXh9Ljd3g+q3Yy2K9aWt+6V1\nze27GgHgsT1vYO/hk4gPp7Aosgtb4p/D4qcuxjXPLMD6335NWiK5JWbkqV7sLj4AxbmRFqxCKtqU\n91RcNOI+fFJrJ1Bsn9RahlfsFcSXymXBKqSe+kKezzEuGrHq7I1o7+33twJxc+sse1XLtSNbnalW\nkYMeW3kTN9eLHbdtfndvv3QF3ZNuzwWBP9H4M6xtfiIvE3TOyHzEnnwFV4/8NE8WOfDeG3D+Wz/N\nBa03pDqxPRs4tbqPdr523FeDEZna5tnoR7Ta/6kqHVp/TzLfvQoB4NE9bziCxc2JY7ihfwMu+NBa\nLP35hb5W3kXXcpndiQbAoYpp1xQM1HPJATbsWSrhi/OVYDG7E3f3HMRt6UdtafEfwr/7lZS5uHV0\nroNqizs+ZkhVFLoBNK8bzvQHR4nyVl728ckaSlhpa4mhfeHtaJ779bznOwC0HfkhZr38cC4rtg0n\ngIGncjuXqQA2Zf+z0t3bj/iws8GIjMWRXegyHsmrAT+VTmB942b8ye/OwGVzrwWgVq6o6rn4mRhV\nyILFDSNDuOxX92P3Sn85GIEY1tmdmZr2AKYC6PLx/vVc150NOyrri/OjQf7umcvxHTgbM/hegSh8\nl/HYFK3roFqJNRkR/0XALLgVrTLT43XG53Y9ZG4aK5f96n7A3gbQI0j9Us+DuGzfBuzDCQw0emvi\nVb7+GIbR9vIGXPCzqdJG48CoT1uGzu/AiJC0K5WJssRxAbGcShvWIAuO1RrsY0ft+OL8NtFQIpFj\npqJNuOudG7Sug5vLRaeOtgqvWvK635PqenhJIgH4D1If2IpZL38NbZSpXT81ktG5q+rHezX9mCLe\nzvmjH93zBoasna1o9PPK/NT2z21tlrKr8Q58dtyL2HjTB3HrFdMdPnPzsbLWvCLG40YhvQGCxKzr\n/tlxL2auQ9Mt2DduaUES4VqjKMNORDcR0UEiShORZ5prtVIrvrjAbpTZnXlB0njsfKxM3obukfnS\nw83rYAbmVOu91paYVCOvi1eDBd3v6cqLJkuPUz2fh8KAHU1PlHfoeW5tXjEzIF/nbn4WswGGV9MP\nu9rGeq1Nl7m5U3mp58FM+YquFuDeWbjvA/+V+32YOwOzWYo17+Hujktw781zcmNriRlozmbFbkh1\nIi6cweKX3v+FvOd0OhfpNswoJaoOXTjgPw+hlijWFfMqgBsAPBjAWCpGpbeMusgkfWMaCpybLXLM\nq9c9j/5hd5WGV8lYr2Cmdnati2tK93va+dpx6eutz1vHZK1/8plxztrppgpG6vpx0VrLauq46en9\nFDq7euSnmPXyw8i5jU4fwWWvrMYjl63B0p9fiBVxyc7A4lIyr7P9e1UVgNv38wuxe/HotdN1XRbb\nG6BoypD/UY0UtWIXQvxCCFFd/ooCqPSW0S9DydHGxoOJZNGZiG47E/M6uAXm3FZiqszJr3a/4rtX\npe735LWyt4/pVDyZkyF+58zlWJm8DfHY+UjD2YbO4fpRrPCPYaJDu3/PDZdgc+OtjhWxEMD/pMdg\nCI24z3hAqxXgioatjp0CkolskPMqTI0oEqhsE5Hse+1Jt6N9eBPed+4xtA9vQk+63bM+UCGuS/v3\ncOk7z+Ky7g9DZHcggayqS1S3qdphHzuqY8uoSyniAaqdiTUDUmUsCXB1uajG+9ieN3zrm3W/J69Y\nhJd65PHhD+Fq8QDePzRq2KzkXQtJvCKBRnyLPo1lW/ryJq2OuW3o+uoaNN/4D8D4aRAgvInJ+OeR\nj6IhInAendHy0wMaQU6VT9z2vK67Uac+kF/XpfV7MF1HbXQCZM2tKNa4a16HsOFp2InoJ0T0quS/\n6/28EREtIaK9RLT3+HH5VrmSFOMbLieliAeoVsJ/1znaSafQwK1qXHY/ve7kZH5P9948BwAcxhPw\nXtnrXKv+rGtARt7zszvx0iVrco2mBzAJd6U+h+8NXZGbtJZt6cMM685kdiew7FVQ1yCmdP0Sfzrx\nUK7qpIlXPZpj8AhyatYr0nE36tYH8uu6tH4P0qCy6TIpBs3rEDY8fexCiI8G8UZCiIcAPARk6rEH\ncc56pNB4gJsvU0dPX6h0TLeKIKA/OXn5eL0+j0pvbyVKpPWZu3v7cedL70Ui+U3lucwfu9IX7eKn\nt9ZyIcooj1pbYhj4wAq02ZqR5BkszXpFss9o1tM3Sx8E9Vuw/wZbmo1crChImWUemtchbLCOvcYo\n5Kbq7u3Hrh88gC34PlrHnMBAfBLu+8EnYS196hXkKrRbjWy8qroquis+nYxG985B3olEI0Lk3idC\noyUI7MFqv0lB0sxLRV7BUPMUTIk1YWAwgbFjGmzX+ypgxgR3g6VRr8jtezUNsb0zViG/BdlkbEQI\nRpSQHBHqoHIQLhON6xA2iuqgREQfB3A/gMkABgH0CSEWer0u9B2USoxfJUHX3auxIvmAtIdp11fX\nlH28qsQbe7ce1We7YOXT0omBAPzaIwFJp3OQeS7VnWGtB+9VB171+rxx2lspIpNXsDJ5Gx63dHgq\nZ0eioPuJqq57S8zA2DENmPfOs1jX+HB+QNiIadUqqifK0kFJCPEDAD8o5hyMf/x2zLlt+FE0R5w+\n3NuGHwWwJjDJmXme/sFEzoVg1lKxt3eb997z0NVzMOcSaTIi2Hv4ZJ7BV7kuipGnFmvUgdGiWTrH\nynCMU+IuuPvsjXh8OD/L2KvOSpDSwaLrvNhQudlOJ5LZBt9XAQfm1p3LpFSwK6YOaFVI31ojbwdW\nTkFVgMrtfOdSo7LNU/Gko1MSIDcmxfh4VYbYzLzUjQkIAN/7jyO+jTpB0Sza5i74rmInoDKQQZfF\nCDJI79Ylyh6IZkMeDCx3rAOGFKV9h2JTPOWTOhmGgF5TDK/jVUayfzCR976FylPdioM1RAn33jzH\ntfywHT/VE03sgVTV9fSrPAlaBhuU8sWrS1S15orUOmzY64DmP1rrqGudijah+Y/Wuq7M/LRl81rJ\n2f9ufWyvaSLTb9vfV0ueemAr4usvQrqrBUdXvR+XvvOscnzJEZEzgssXzpTWH7djlgkoFLNZtAy/\nSXNBy2CLTdozFwRLt/R5doligocNez0wuxMN19+f10Cj4fr7gdmdriszP6tAr5Wc/e/mY1lNE1ly\nju/Vp9kI2dIj1ivpxzSCHXPbcIukUJaVmBHFp35vmtT4TWg2FK9yMphIBlJnJbACcRasCqAJzYa2\nIbYuCFSYXaKY0sCGvV7IJsWgazCvmYbbyszPKtCtObJspWce76dZtK/Vp6IRslvST4QoZ2RlhbIm\nNBt5Rvbujkukxnf1oot9NYpWTVh+kuaCLIthGmar1t9axsILHQlotdVhKjW6Ls2g4OBpiChEFeGm\nSVb1MZXdlPbep3ZVjKw4FAC0PqUI7MoaS/sxBi5JPypGhHAkOulcP9UxX9q6X8sPH0QV0ULzDGQU\nq4jx+jz15luvRL8HNuwhoZgfj8o4+VWf6MowrRPQFU2TMAXOEhP28rW+jYFLI2RzspEZ3mIkfVbM\n17tVxDSRTViFTtJBGIpi/fVuyiLVRB9mgpaO6hA+w35ga11qYUvx4wlyFWhin4C+MXyTo6OQWb42\nSoS0ENrvazWGnxl3I+6i/4NGMZrwYjZCvvKiydi445ByNR1UHX779RsfM3B2OJXXTFo2YfmdpIMu\nfVtsGWvVgqDWg6WFXudK9HsIl2G3Z/CZFeKA0Bv3Uv14gloFmtgnIFX972ejH8Hf+TAEdmP4nTOX\n40xjCl8Zsw3jk2/lGiE3zfqEI+vVjpsB83tz26+fzutVk/Sa7QcdrwUQ+Da/2JZypVgQVJpidsSV\n6PdQVEmBQilZSYF7Z0m33xg/LRMwDDGqlG1rs4dqQFUOAMiM1a8hsGa7qs5p/fxeJQXcVpZBp9mr\ncLtG9rGOaYhIC5oV+71Xorl7NVPM/RXk76YsJQWqjjotqg/UTuNer4bVfvDq6gS46+dlY3AzYOXy\nlepmvyaSI8rPXm07tVqnmB1xJXYw4TLsioBZ2IvqA4X9eApalRUQw7C3oTMihGRa7mf2MyYdWV2L\nTVNezMRSLl+pbJL2S73JCUtNse6Uck+U4dKx12lRfRM/umc/WaU5zBjG6SOAZpcbWRs6q1G3Jr74\nHZOOQbV7GovRe5ciCUhGx9w23HhpG3QSWyc0GzXV1rFWqbX2meEy7LM7M2U+LRmWXPZTTkG1RRSN\ngeM/Vk+cXqtqa+KL3zHpGNTTNv9zoXVmgPLd3N29/djy4hHHpGQnZkSxetHFNdPWsZYp5ndTCcLl\nigG4QpwmBbkVFLGKpvib6O7tl/7IvVbVVsOt8iurzqHjslAlUxVyQ5bLV7pxx6G8XY0Me62VajQw\nYQvA1lLcIXyGndGiIJ+hS9KPKoCoEwi0SsdkCMBZ5hbObFc7pVhNl+Pm1nExVXutlUpkWzKjhMsV\nw2hTkFthwSrERWPeU2Yykduq2stVHCXyDBSq/O1mXOG/112H+yy1Xap9q+yGjoup2oOjQZcRZvzB\nK/Y6pSC3wuxObOg5iNuGH81LJupJtytrmHfMbcPSLX3KU8aMqLb6w2oYdDs1VRs67onlC2di+bb9\nSndMNQftTCqRbcmMwoa9jinErTDnuiW4+snfR2JYXy/fpnDHmH5itwQjO/2DCSx/fH8uLV+nU1O1\noOueMP/f2jrQbKhdK7VWKpFtyYzChr2MhCGYVMhKX6d2iP3vbr1ErbVWrJS6sFKx+ElwqqVAnYxa\nSZgLK2zYy0SYgkl+jY7XZCD7+5UXTfas6SKjmrf69eSeCGO9mFqiKMNORBsBLAIwDOBXAP63EGIw\niIGFjUqU7qwmvCYD2d/nvfc8V/+8DHumadAUs+uqN/dEre86apliVTHPApglhJgN4D8B3Fn8kMJJ\nPa3WgqJjbpt2Y2mTUta0Kyhb10KtZS8ytUtRhl0I8YwQIpV9uAdA+IuyFEi50tHDhswYGhG1gNKe\naRokxUr4ai17kaldgvSx/xmALao/EtESAEsAYPr06QG+bW3AwaTCUPlq/bTtC4ogdl3snmDKgadh\nJ6KfAJgi+dNXhBBPZY/5CoAUgMdU5xFCPATgISBTj72g0dYwHEwqHJUxLPdEWW8+ci/CoPIKK56G\nXQjxUbe/E9FnAfwxgAWiEl07aojAV2t12gYQqMxEybuuUcKk8gojxapirgWwAsBHhBDxYIbEaFHH\nbQBNvCbKoFeUvOsapd5VXtVOsT72bwEYA+BZyhSP3iOE+IuiR8V4oyihi+fWKg17PW2dS7WiZB95\nBlZ5VTdFGXYhxP8KaiCMT3y2Aay3rTOvKEsLxxuqG67uWKuo2v0pnq+3anu8oiwtrMmvbtiw1yo+\n2wDWm6HjvIHSwpr86oZrxdQqph9dUxVTb1vnwBUsdgXShdcA//VMVSuSSh1TcYs31FM8pxqhSigU\n582bJ/bu3Vv2961n7D52wFlhMWzYjcuVF03GzteO+zc2dgWSDCNWVf11K/l91+NvrVwQ0T4hxDzP\n49iw1w/1vIoqytjcO0vaEtDB+GnAsleLHGkwzF/3vHSH1tYSK6oZic5vqFTvzegbdnbF1BH1LNUr\nSiWjUiAVelwZKEVMRVdZVW/xnGqEg6dMXVCUsVEpkAo9rgyUInisq6ziwHXlYcPO1AVFGRuZAsmO\niyKpEpRCjqg7ObIUsvKwYWfqgqKMzezOTGB0/DQAlPl33p/nP66iwClQGjmi7uTIUsjKw8FTpm6o\n5+BxELDapfJw8JRhbNRz8DgIuAha7cCGnWEYbXhyrA3YsDMMU/Owmy0fNuwMw9Q09Va5VAdWxTAM\nU9PUW+VSHdiwMwxT03CmqxM27AzD1DSc6eqEDTsTOrp7+zF/3fO4YOXTmL/ueXT39ld6SEwJ4UxX\nJxw8ZUIFB9LqD9bXO2HDzoQK7nVan7C+Pp+iXDFE9LdEdICI+ojoGSJqDWpgDFMIHEhjmOJ97BuF\nELOFEHMA/BBA9ZS3Y+oSDqQxTJGGXQjxjuXhWADlryjGMBY4kMYwAfjYiejrAP4UwGkAVxY9IoYp\nAg6kMYxG2V4i+gmAKZI/fUUI8ZTluDsBNAkhVivOswTAEgCYPn36pYcPHy540AzDMPVI2ZtZE9F0\nAD8SQszyOpbrsTMMw/hH17AXq4q50PLwegCvFXM+hmEYpniK9bGvI6KZANIADgP4i+KHxDAMwxRD\nUYZdCHFjUANhGIZhgoFrxTAMw4QMNuwMwzAhgw07wzBMyGDDzjAMEzLYsDMMw4QMNuwMwzAhgw07\nE24ObAXunQV0tWT+PbC10iNimJLDjTaY8HJgK7D9DiCZrcV++kjmMQDM7qzcuBimxPCKnQkvz60d\nNeomyUTmeYYJMWzYmfBy+qi/5xkmJLBhZ8LL+Kn+nmeYkMCGnQkvC1YBhq0lnhHLPM8wIYYNOxNe\nZncCizYB46cBoMy/izZx4JQJPayKYcLN7E425EzdwSt2hmGYkMGGnWEYJmSwYWcYhgkZbNgZhmFC\nBht2hmGYkMGGnWEYJmSwYWcYhgkZbNgZhmFCBgkhyv+mRMcBHC77GzuZBOBEpQdRALU6bqB2x16r\n4wZqd+y1Om6gdGN/rxBistdBFTHs1QIR7RVCzKv0OPxSq+MGanfstTpuoHbHXqvjBio/dnbFMAzD\nhGZmfn0AAANkSURBVAw27AzDMCGj3g37Q5UeQIHU6riB2h17rY4bqN2x1+q4gQqPva597AzDMGGk\n3lfsDMMwoaOuDTsR/S0RHSCiPiJ6hohaKz0mXYhoIxG9lh3/D4iopdJj0oGIbiKig0SUJqKaUDwQ\n0bVEdIiIfklEKys9Hl2I6NtE9BYRvVrpsfiBiKYR0U4i+nn2t/LFSo9JByJqIqIXiWh/dtxrKjaW\nenbFENG7hRDvZP//DgAfEEL8RYWHpQURXQPgeSFEiojWA4AQ4ssVHpYnRPQ7ANIAHgTwN0KIvRUe\nkitEFAXwnwCuBnAUwEsAPiWE+HlFB6YBEX0YwBkAjwghZlV6PLoQ0fkAzhdCvExE7wKwD0BHtV9z\nIiIAY4UQZ4jIALALwBeFEHvKPZa6XrGbRj3LWAA1M8sJIZ4RQqSyD/cAqIkOzUKIXwghDlV6HD64\nHMAvhRCvCyGGAXwfwPUVHpMWQogXAJys9Dj8IoQ4JoR4Ofv//wPgFwDaKjsqb0SGM9mHRva/itiU\nujbsAEBEXyeiIwBuAVCrXY7/DMCPKz2IkNIG4Ijl8VHUgJEJC0Q0A8BcAP9R2ZHoQURRIuoD8BaA\nZ4UQFRl36A07Ef2EiF6V/Hc9AAghviKEmAbgMQB/XdnR5uM19uwxXwGQQmb8VYHOuBnGCyIaB+AJ\nAEttu+uqRQgxIoSYg8wO+nIiqogLLPTNrIUQH9U89DEAPwKwuoTD8YXX2InoswD+GMACUUXBEh/X\nvBboBzDN8nhq9jmmhGR91E8AeEwI8WSlx+MXIcQgEe0EcC2AsgevQ79id4OILrQ8vB7Aa5Uai1+I\n6FoAKwAsFkLEKz2eEPMSgAuJ6AIiagTwSQA9FR5TqMkGIR8G8AshxN9Xejy6ENFkU51GRDFkAu4V\nsSn1rop5AsBMZFQahwH8hRCiJlZjRPRLAGMAvJ19ak8tKHqI6OMA7gcwGcAggD4hxMLKjsodIvoY\ngPsARAF8Wwjx9QoPSQsi+h6AP0Sm0uBvAKwWQjxc0UFpQETtAP4NwCvI3JsAcJcQ4keVG5U3RDQb\nwHeR+Z1EAGwVQqytyFjq2bAzDMOEkbp2xTAMw4QRNuwMwzAhgw07wzBMyGDDzjAMEzLYsDMMw4QM\nNuwMwzAhgw07wzBMyGDDzjAMEzL+P/pYrC0m3x6eAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x118045a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "An MMD at least this large would occur with probability 0.13 under the hypothesis that the two distributions are the same.\n"
     ]
    }
   ],
   "source": [
    "X = np.random.randn(500,2)\n",
    "Y = np.random.randn(100,2) + 0.\n",
    "\n",
    "plt.scatter(X[:,0], X[:,1])\n",
    "plt.scatter(Y[:,0], Y[:,1])\n",
    "plt.show()\n",
    "\n",
    "mmd, p_value =  MMD_test(X, Y, gaussian_kernel, 1)\n",
    "print \"An MMD at least this large would occur with probability \" + str(p_value) + \\\n",
    "\" under the hypothesis that the two distributions are the same.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hilbert Schmidt Independence Criterion: are X and Y independent?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use a very similar methodology for independence testing.\n",
    "\n",
    "Suppose we are given samples $(X_i,Y_i)$, $i=1,\\ldots,N$ draw from the joint distribution $\\mathbb{P}_{XY}$. We want to know whether $X$ and $Y$ are independent. Equivalently, we can ask whether $\\mathbb{P}_{XY} = \\mathbb{P}_{X}\\mathbb{P}_{Y}$.\n",
    "\n",
    "It can be shown that if we have kernels $k$ and $l$ over $\\mathcal{X}$ and $\\mathcal{Y}$ respectively, then $k \\otimes l$ is a kernel on $\\mathcal{X}\\times \\mathcal{Y}$. We can write the mean embedding of $\\mathbb{P}_{XY}$ with this kernel as $\\mu_{XY}$. \n",
    "\n",
    "It can also be shown that if $\\mathbb{P}_{XY} = \\mathbb{P}_{X}\\mathbb{P}_{Y}$ then $\\mu_{XY} = \\mu_{X}\\otimes \\mu_{Y}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have empirical estimates for the above quantities:\n",
    "\n",
    "$$\n",
    "\\tilde{\\mu}_{XY} = \\frac{1}{N} \\sum_{n=1}^N k(X_n, \\cdot) l(Y_n, \\cdot)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\tilde{\\mu}_{X}\\otimes\\tilde{\\mu}_Y = \\frac{1}{N} \\sum_{n=1}^N k(X_n, \\cdot)  \\frac{1}{N} \\sum_{m=1}^N l(Y_{m}, \\cdot)\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our test statistic will be \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\| \\tilde{\\mu}_{XY} - \\tilde{\\mu}_{X}\\otimes\\tilde{\\mu}_Y \\|^2 &= \\langle \\tilde{\\mu}_{XY}, \\tilde{\\mu}_{XY}\\rangle - 2 \\langle \\tilde{\\mu}_{XY}, \\tilde{\\mu}_{X}\\otimes\\tilde{\\mu}_Y\\rangle + \\langle \\tilde{\\mu}_{X}\\otimes\\tilde{\\mu}_Y, \\tilde{\\mu}_{X}\\otimes\\tilde{\\mu}_Y\\rangle \\\\\n",
    "&= \\frac{1}{N^2} \\sum_{n,m} k(X_n,X_m)l(Y_n,Y_m) - \\frac{2}{N^3}\\sum_{n,m,s} k(X_n,X_m)l(X_n,X_s) + \\frac{1}{N^4}\\sum_{n,m,s,t} k(X_n,X_m)l(X_s,X_t)\\\\\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll write $K$ and $L$ for the matrices with entries $K_{ij} = k(X_i, X_j)$ and $L_{ij} = l(Y_i, Y_j)$.\n",
    "\n",
    "We will also write $M_{++}$ for the sum of all of the entries of the matrix $M$, $\\circ$ for the Hadamard (aka element-wise) product (e.g. $K\\circ L$) and $M_{+}$ for the vector with entries $(M_{+})_i = \\sum_n M_{in}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this abbreviated notation, we can express the above quantities more succinctly:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\| \\tilde{\\mu}_{XY} - \\tilde{\\mu}_{X}\\otimes\\tilde{\\mu}_Y \\|^2 &= \\frac{1}{N^2} (K\\circ L)_{++} - \\frac{2}{N^3} (K_{+}\\circ L_{+})_{+}  + \\frac{1}{N^4} K_{++} L_{++}\\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "This test statistic is known as HSIC (Hilbert Schmidt Independence Criterion)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we won't expect in the finite-sample case that this statistic will actually be equal to $0$ even if $X$ and $Y$ are independent. We need to see if it is _significantly far away from zero_ under the hypothesis of independence. We can simulate draws under the hypothesis of independence by applying a permutation to the indices of one of the $X_i$ or $Y_i$.\n",
    "\n",
    "That is, we can generate lots of permutations $\\pi$ of the integers $1,\\ldots,N$ and calculate HSIC using the dataset $(X_{\\pi(i)}, Y_i)$ and get a p-value for our actual sample of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "K = np.random.randn(5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17.860681577107975"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "( K.sum(axis=1) * K.sum(axis=1) ).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 362,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def HSIC(X, Y, kernel_X, parameters_X, kernel_Y, parameters_Y, null=False):\n",
    "    if null is True:\n",
    "        X = np.random.permutation(X)\n",
    "    K = kernel_X(X, parameters_X)\n",
    "    L = kernel_Y(Y, parameters_Y)\n",
    "    N = len(X)\n",
    "    \n",
    "    hsic = ( (K*L).sum() / float(N**2) ) - ( 2 * (K.sum(axis=1) * L.sum(axis=1)).sum() / float(N**3) ) + ( K.sum() * L.sum() / float(N**4) )\n",
    "    return hsic\n",
    "\n",
    "def HSIC_test(X, Y, kernel_X, parameters_X, kernel_Y, parameters_Y):\n",
    "    hsic = HSIC(X, Y, kernel_X, parameters_X, kernel_Y, parameters_Y, null=False)\n",
    "    \n",
    "    n_samples = 1000\n",
    "    null_distribution = [False for _ in range(n_samples)]\n",
    "    for i in range(n_samples):\n",
    "        null_distribution[i] = HSIC(X, Y, kernel_X, parameters_X, kernel_Y, parameters_Y, null=True)\n",
    "    null_distribution = np.array(null_distribution)    \n",
    "    p_value = ( null_distribution > np.array(hsic) ).sum() / float(n_samples)\n",
    "    return p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 368,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF05JREFUeJzt3X+IZWd9x/HPN+NoxygdS4aaTDJuaO3aNBtcOkQhhdao\n3fRn1i2CwbaIha1QQUFWNi40ShsSWJCKCu3SiBWDGptkTTVlE8mCNXTTzGY3JnF3JVgSM0pda8Zf\nOzSzm2//mJnNZHLPvefc85xznuc57xcEMjN373nO3Dvf+5zv832+x9xdAIB8XND1AAAAYRHYASAz\nBHYAyAyBHQAyQ2AHgMwQ2AEgMwR2AMgMgR0AMlM7sJvZL5nZf5nZo2b2hJl9LMTAAADjsbo7T83M\nJF3o7j83s0lJ35T0AXc/UvRvLrroIt+yZUut4wJA3xw9evRH7j4z6nEvq3sgX/1k+Pnal5Nr/w39\ntNiyZYsWFhbqHhoAesXMnirzuCA5djObMLPjkn4o6X53fyjE8wIAqgsS2N39nLu/UdKlkq42sys3\nP8bMdpvZgpktnD59OsRhAQADBK2KcfclSYclXTfgZwfcfd7d52dmRqaIAABjClEVM2Nm02v/PyXp\n7ZJO1n1eAMB4ai+eSrpY0r+Y2YRWPyjucPevBnheAMAYQlTFfEvS9gBjAQAEEGLGDmTv4LFF7T90\nSt9fWtYl01Pas2Ordm6f7XpYwEAEdmCEg8cWdeNdj2l55ZwkaXFpWTfe9ZgkEdwRJXrFACPsP3Tq\nfFBft7xyTvsPnepoRMBwBHZghO8vLVf6PtA1AjswwiXTU5W+D3SNwA6MsGfHVk1NTrzoe1OTE9qz\nY2tHIwKGY/EUGGF9gZSqGKSCwA6UsHP7LIEcySAVAwCZIbADQGYI7ACQGXLsQMRoZYBxENiBSNHK\nAOMiFQNEilYGGBeBHYgUrQwwLgI7EClaGWBcBHYgUrQywLhYPAUiRSsDjIvADkSMVgYYB6kYAMgM\ngR0AMkNgB4DMENgBIDMsniJKZXqk0EcFGIzAjuiU6ZFCHxWgGKkYRKdMjxT6qADFCOyITpkeKfRR\nAYoR2BGdMj1S6KMCFCOwIzpleqTQRwUoVnvx1Mwuk/Q5Sb8qySUdcPdP1H1e9FeZHin0UVlFZRAG\nMXev9wRmF0u62N0fMbNXSzoqaae7f7vo38zPz/vCwkKt4wJ9t7kySFq9arll1zaCe6bM7Ki7z496\nXO1UjLv/wN0fWfv/n0k6IYl3FUo5eGxR19z6gC7f+zVdc+sDOnhsseshJYPKIBQJWsduZlskbZf0\n0ICf7Za0W5Lm5uZCHhaJoha9HiqDUCTY4qmZvUrSnZI+6O4/3fxzdz/g7vPuPj8zMxPqsEgYM856\nqAxCkSCB3cwmtRrUb3f3u0I8J/LHjLMeKoNQpHZgNzOTdJukE+7+8fpDQl8w46xn5/ZZ3bJrm2an\np2SSZqenWDiFpDA59msk/YWkx8zs+Nr3PuLu9wZ4bmRsz46tA6s6mHGWl8odlijLbFftwO7u35Rk\nAcaCnqEWvR9YJG8f3R3RqVRmnBjfsEVyXvtm0FIAQKNYJG8fgR1Ao1gkbx+BPVLsyEQuKMtsHzn2\nCLHYhJywSN4+AnuEWGxCKLGUGbJI3i4Ce4RYbEIIXPn1Fzn2CLHYhBDoxdNfBPYIsdiEELjyqy/V\nIgZSMRFisQkhXDI9pcUBQZwrv3JSTmUR2CPFYhPqohdPPSkXMRDYgUxx5VdPyqksAjuQMa78xpdy\nKovFU2BNqgtlaEbKRQzM2AGlvVCGZqScyiKwA0p7oQzNSTWVRSoGUNoLZcBmBHZA7PZFXkjFoLc2\nNsiafuWkJi8wrTzv539eZaEslmZbgERgR09tXix99syKJidM01OT+snySqXgzMIrYkNgjxCzv+YN\nWixdOee68BUv0/Gbfr/2c7Hwii4R2CPD7K8dIRdLWXhFbFg8jQytVtsRcrF01HOx8QltY8YeGWZ/\n49mcvnrLG2Z0+OTpwnRWyAZZw56LKzB0gRl7ZCi7q249eC4uLcu1Gjw/f+TpF319412PvWimvHP7\nrG7ZtU2z01MySbPTU7pl17axgu2w5+IKDF1gxh4ZWq1WNyh4bjZoMTPkrsKi5ypzBcZiOUIjsEcm\n5f4UXSmbpuoinTWqQyCpGjSBwB6hVPtTdKUoeA56XNtGXYFRKlmMK5nxBcmxm9lnzOyHZvZ4iOcD\nqhjUXnWzrtJZo3L5LJYPNmjdZPM6CYqFmrF/VtKnJH0u0PMBpQ1KX42qiml7fEXHTvlmDk3iSqae\nIIHd3b9hZltCPBcwStEleop/8CyWD8aVTD3k2NGqunnT3BYbWSwfjCuZeloL7Ga2W9JuSZqbm2vr\nsIhIiKCc4yV6qlcbTeJKpp7WNii5+wF3n3f3+ZmZmbYOi4iE2KzDJXr7umiJEHIDWR+RikFrQgTl\nokv0C8x08Ngif/iBdZn64kpmfKHKHb8g6T8lbTWzZ8zsr0I8L/ISol1CUWnjOXfK4RpAS4Q0BQns\n7n6Du1/s7pPufqm73xbieZGXQUG5at50/RJ9wuwlPyPghEfqK000AUNrQuVNd26f1fPuA39GwAmL\npnRpIseOVoXKm6ZcDpfSVnmqU9LEjB1JCpHW6ULorfJNV6xQnZImZuxIUiwbe6re4CNkHX5bFStU\np6SHwI5klQ04TaU+BgXWzx95+vzPBwXakIuROW7WQhikYpC1QamPD37puN74sftqpy2q3OBjXcjF\nyNAVK9ybNR8EdmStKPguLa/Urnsf5wYfIdcGQn5I0CY3LwR2ZG1Y8F1eOacP3fHo2MGrbADd+LiQ\ni5EhPyTYiJQXAjuyNir41tmx2vUNPup+SGxMvRTdgYp9AWli8RRZG1SHvdm4C47j3OAjdCXLuBUr\nm8dRJIV9AXgpAjuyth70PvZvT+jZMyuFjxt3Zro5sB48tqjDJ08XPj50Jcu4FT9lFn5T2BeAwQjs\nyN568D14bFEfuuNRnRvQjiDEzLTMbDxkJUud2f+w45kU/Y5YDEdgx0gpbYEfZn3MTW2RLzMbD9kK\noc7sv2gcs9NTenDvtQP/TS7vgz5g8RRD5VYG1+QW+TKz8aIF11/839nKv9M6s/+qFTW5vQ9yx4w9\nc3VnWTnubmxqi3yZ2fj6cffd/Zh+8dwLv9f1uvqNjwlxvCJVWzLk+D7IGYE9YyEqMOjHXV6VTohn\nnnvpwmXVQFm382KVDzjeB2khFZOxEJtO6MddXtk0z/5DpzS4m3y1QNlm50XeB2lhxp6xELOslPpx\nl0k7Nb0AWGYWPOz3XzVQDjteyHNN6X0AAnvWQlRgxNIed5Qyaacub8y8UdHrYlKwQNnERigp/vcB\nVpkX3GKsSfPz876wsND6cftm0O7CqcmJLG+UcM2tD4ws3yvzmDYMel1M0rvfPKe/37ktyDFiOVes\nCnX1ZGZH3X1+1OOYsWesT7OsMmmnWBYA23hdYjlXdHOlSGDPXF/uflMm7RTTfVKbfl1iOte+66JU\nlKoYZKHMhps6bW5TuwlFm/eETe1307Yurp6YsSMLZdIb46ZAYll03TieUefQVhoutt9NjLq4emLx\nFElro39JTAuRsS2Ix/S7iVXI16zs4impGCSrrf4lMS1Exnano5h+N7FqcyPZOlIxSFZbi1IxLUTG\nFkhj+t2si7ELZdtFDMzYkay2glzohcg6i411t/aHXuhsc5G2DLpQriKwI1lt9S8JeSldN/DUrewJ\nHfS6SDMME1uqqitBUjFmdp2kT0iakPTP7n5riOcFhmmzf0moS+m66aM61S5Npa5i2isRW6qqK7UD\nu5lNSPq0pLdLekbSw2Z2j7t/u+5zA8OkuLM2ROAZN5D2IejFmPPvQogZ+9WSnnT370qSmX1R0vWS\nCOxoXEyzxTKLdl0Gnj4EPbpQrgoR2GclfW/D189IetPmB5nZbkm7JWlubi7AYYFmVamuKLtRZ5zA\nE6rKI+agF+ocU7yKa0Jr5Y7ufkDSAWl1g1Jbx8VgMZaExaTqjsqy+euqgSfkzs66Qa+p90wTLYb7\n/l4OEdgXJV224etL176HSLENfLSqC41V8tdVAk/oBc9xg16T7xnupxpeiHLHhyW93swuN7OXS3qX\npHsCPC8aMk5JWN8aPVVdaGyq9DKWBc8mywhjOcec1A7s7n5W0vslHZJ0QtId7v5E3edFc6r+IfVx\n00fVQN3URp1Y7jXaZPCN5RxzEmSDkrvf6+6/4e6/5u43h3hONKfqH1IfN31UDdRNbdSJZWdnk8E3\nlnPMCb1ieqhqdURXl8pdLvCOs9DYxKJdLFUeTVbUxHKOOaFtb09VCZpdtGaNrT1tykJ9QFJJ1b2y\nbXsJ7BipiyBLn+8w+IDMCzezRjBdXCrnXinR1uyXUsJ+IrCjlLY3feS8/b3NfQS5f0BiMNr2Iko5\nV0qMqjIKuWeAUsJ+IrAjSrH1+Q5p2Cw69J6BnD8gUYxUTAL6Wo2Qa8+PYWmmJloISJQS9g0z9sj1\ncddn7obNopvIie/cPqs9O7bqkukpfX9pWfsPneL9kzkCe+Ry3PXZt74zm89XUmGaqYmcOJOD/iEV\nE7ncqhr61lmy6Hxv2bVtYD1+Ezs8KXnsHwJ75HIr+ysKMh+954lW88Cb1y3e8oYZHT55OvjxqwbV\nJnLiuU0OMBqBPXKx3PUm1AJuUTBZWl7R0vKKpOZn8YNm0Z8/8vT5n4c8/jhBNfSicVeTg74u+seA\nHHvkYij7C5mjLRtMmlxHGDSLbur4MdSRd1HySF6/W8zYE9B12V+IHO367G1xaVkmqUyHoqZSBWWf\nN8TxY7ji6qLkkbx+twjsGGmcG3Nszl/feXTx/B+6S+eD++z0lM48d1bPnll5yfM0NastSk00cfzY\n68ibSpeQ1+8WgR0jVcnRDspf337k6ZfM0NeD+oN7ry3sQFh3VlsUtAbNojcLOavu+oqrqDJn4akf\nv+gDN+TaQm6L/qkhx46RquRoB12CF6Vd1mdvTawjDMvxDjren795Lsv2BVJxWuQLD32vsT0SXbcy\n6Nteic2YsWOkKumEKpfaG2dvoWe1o3K8Xc+i21T0mpwruBdDiHRJlymovu2VGITAHqEYy8TKBsKi\nS/DNC6ZNz97I8b6g6DWZMBsY3EOlS7r68GThllRMdFIvEyu6BH93y6mOGMoMY1H0mtzwpsuy7PzI\nhzoz9uikPtuIpQokhjLDWAx7TeZf9yudv1ahsXDLPU+jc/nerw1cbDRJ/33rH7U9nKTFmNJC83K+\nzyv3PE0Us41w+rRAihfEctXYJQJ7ZEghAPX1/UOdwB4ZZhsA6iKwR6jvsw0A9VDuCACZqTVjN7N3\nSvqopN+UdLW7U+rSkL5WeMR63rGOC5Dqp2Iel7RL0j8FGAsK9HWLdKznHeu4JD5wsKpWKsbdT7h7\nundVTkRMN7Rus7lSTOe9UazjSn3XMsIhx56AWLZItx04YjnvssfvelyxfuCgfSMDu5l93cweH/Df\n9VUOZGa7zWzBzBZOnz49/oh7KJa+J20HjljOu+zxi77f1lVOrB84aN/IwO7ub3P3Kwf895UqB3L3\nA+4+7+7zMzMz44+4h7rubb2u7cARy3lvVmVcbV7lxPpBiPZRx56AWDYttd3uIJbzrjOuJpu6jboF\noRTHByHaV6sJmJm9Q9InJc1IWpJ03N13jPp3NAFLU87NlZrSVFO3otfiz357VodPno7qgxDhtNIE\nzN3vlnR3nedAOmKdQceszlXOsNLFoiuBwydP68G919Ye98Fji/roPU9oaXn1JuOveeWkbvqT3+K1\nTgSpGFRCu4Nqxm3qNqpWvsn1joPHFrXny49q5fkXrjWePbOiPf/66PnjI26UOyJLsdzMeNwbdY+q\nQGpyoXT/oVMvCurrVs45pZOJYMaO7MS2M3Scq5xRM/Im2zsPm/VTOpkGZuzITg4bdUbNyMe9Eqhz\n7FE/QzyYsSNJwxYWc9ioU2ZG3tR6x54dW1+SY5ekyQmjdDIRBHZEoUrzqlGplhxuL9hlBdL6MaiK\nSRc3s0bnqtbHX3PrAwMD9+z0lB7cey319sgWN7PuqRTbtlbdnTkq1UK9PfqOwJ6R2KpByqqaEy+T\naqHeHn1GVUxGUq0GqVqTHWtzMCAWBPaMpFoNUjVQN1nqB+SAVExGUq0GGScnTqqlvBTXXVAPgT0j\nTe5GbBqBuhmprrugHgJ7RqgGwWZN9oOviiuH9hDYE1X0R8LMFxvFsu7ClUO7WDxNEHejR1mx3C4v\n1YqtVBHYE8QfCcqKpTQ0liuHviCwJ4g/EpQVS2loLFcOfUGOPUGpljWiGzGsu6RcsZWipGbssdwV\np2uxXF4DZcVy5dAXyczYWVV/AWWNSFEMVw59kUxgj6keNwb8kQAokkwqhgVDACgnmcDOqjoAlJNM\nYGfBEADKSSbHzoIhAJSTTGCXWDAEgDKSCuxAjOhaiNgQ2IEa2F+BGNVaPDWz/WZ20sy+ZWZ3m9l0\nqIEBKaAhG2JUtyrmfklXuvtVkr4j6cb6QwLSwf4KxKhWYHf3+9z97NqXRyRdWn9IQDrYX4EYhaxj\nf6+kfw/4fED02F+BGI1cPDWzr0t67YAf7XP3r6w9Zp+ks5JuH/I8uyXtlqS5ubmxBgvEhv0ViJG5\ne70nMHuPpL+W9FZ3P1Pm38zPz/vCwkKt4wLoDiWe3TCzo+4+P+pxtcodzew6SR+W9LtlgzqAtFHi\nGb+6OfZPSXq1pPvN7LiZ/WOAMQGIGCWe8as1Y3f3Xw81EABpoMQzfsl0dwQQB0o840dgB1AJJZ7x\no1cMgEoo8YwfgR1AZbTQjhupGADIDIEdADJDYAeAzJBjBzLH9v/+IbADmdkYyH95alK/eO6sVs6t\n9oRi+38/kIoBMrLex2VxaVkuaWl55XxQX8f2//wR2IGMDOrjMgjb//NGYAcyUjZgs/0/bwR2ICNl\nAjbb//NHYAcyMqiPy+QFpte8clImaXZ6Srfs2sbCaeaoigEyQh8XSAR2IDv0cQGpGADIDIEdADJD\nYAeAzBDYASAzBHYAyAyBHQAyY+4++lGhD2p2WtJTY/zTiyT9KPBwYtOHc5T6cZ59OEepH+cZyzm+\nzt1nRj2ok8A+LjNbcPf5rsfRpD6co9SP8+zDOUr9OM/UzpFUDABkhsAOAJlJLbAf6HoALejDOUr9\nOM8+nKPUj/NM6hyTyrEDAEZLbcYOABghucBuZn9nZt8ys+Nmdp+ZXdL1mEIzs/1mdnLtPO82s+mu\nx9QEM3unmT1hZs+bWTIVB2WY2XVmdsrMnjSzvV2Ppwlm9hkz+6GZPd71WJpiZpeZ2WEz+/bae/UD\nXY+pjOQCu6T97n6Vu79R0lcl/W3XA2rA/ZKudPerJH1H0o0dj6cpj0vaJekbXQ8kJDObkPRpSX8g\n6QpJN5jZFd2OqhGflXRd14No2FlJH3L3KyS9WdLfpPBaJhfY3f2nG768UFJ2iwTufp+7n1378oik\nS7scT1Pc/YS7n+p6HA24WtKT7v5dd39O0hclXd/xmIJz929I+nHX42iSu//A3R9Z+/+fSTohKfpm\n90neaMPMbpb0l5J+IuktHQ+nae+V9KWuB4FKZiV9b8PXz0h6U0djQSBmtkXSdkkPdTuS0aIM7Gb2\ndUmvHfCjfe7+FXffJ2mfmd0o6f2Sbmp1gAGMOse1x+zT6qXg7W2OLaQy5wnEzsxeJelOSR/clDWI\nUpSB3d3fVvKht0u6VwkG9lHnaGbvkfTHkt7qCdekVngtc7Io6bINX1+69j0kyMwmtRrUb3f3u7oe\nTxnJ5djN7PUbvrxe0smuxtIUM7tO0ocl/am7n+l6PKjsYUmvN7PLzezlkt4l6Z6Ox4QxmJlJuk3S\nCXf/eNfjKSu5DUpmdqekrZKe12qHyPe5e1azITN7UtIrJP3v2reOuPv7OhxSI8zsHZI+KWlG0pKk\n4+6+o9tRhWFmfyjpHyRNSPqMu9/c8ZCCM7MvSPo9rXY+/B9JN7n7bZ0OKjAz+x1J/yHpMa3GHEn6\niLvf292oRksusAMAhksuFQMAGI7ADgCZIbADQGYI7ACQGQI7AGSGwA4AmSGwA0BmCOwAkJn/BxbA\n5BDITPgvAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1176af0d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of HSIC being its size or larger under the hypothesis that X and Y are independent is 0.295\n"
     ]
    }
   ],
   "source": [
    "X = np.random.randn(100,1)[:, None]\n",
    "Y = np.random.randn(100,1)[:, None]\n",
    "\n",
    "p_value = HSIC_test(X, Y, gaussian_kernel, 1, gaussian_kernel, 1)\n",
    "\n",
    "plt.scatter(X[:,0], Y[:,0])\n",
    "plt.show()\n",
    "print \"The probability of HSIC being its size or larger under the hypothesis that X and Y are independent is \" + str(p_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 372,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF7xJREFUeJzt3X9sXWd9x/HPN+4tOIBwWKOxug2pNBREKW2EBUX5Y6PA\nUsZKQwsbTDAxJmVIQwJUhaVUokUw1VM0pklMYtVA/LEICksJhTKFVqmE1K0Mm6Q/0jaIH2qpi0RQ\na2CNoY7z3R+xXfv6nHt+Peeec577fkmVsH197nOBfvyc7/k+z2PuLgBAPDY1PQAAQFgEOwBEhmAH\ngMgQ7AAQGYIdACJDsANAZAh2AIgMwQ4AkSHYASAy5zXxphdccIFv3769ibcGgM6anZ39pbtvzXpd\nI8G+fft2zczMNPHWANBZZvZ4ntdRigGAyBDsABAZgh0AIkOwA0BkCHYAiAzBDgCRaaTdEQBGyeFj\nczpw5KSeml/QhRPj2rd7h/bsnKzt/Qh2AKjR4WNzuvGOh7SwuCRJmptf0I13PCRJtYU7pRgAqNGB\nIydXQ33FwuKSDhw5Wdt7EuwAUKOn5hcKfT8Egh0AanThxHih74dAsANAjfbt3qHx3ti67433xrRv\n947a3pOHpwBQo5UHpHTFAEBE9uycrDXI+1GKAYDIEOwAEBmCHQAiQ7ADQGQIdgCIDMEOAJEh2AEg\nMgQ7AESGYAeAyBDsABAZgh0AIsNeMQCwxrCPsatD5Rm7mV1sZvea2SNmdsLMPhJiYAAwbCvH2M3N\nL8j1/DF2h4/NNT20QkLM2M9IusHdf2BmL5E0a2Z3u/sjAa4NAJUUmYEPOsauS7P2ysHu7j+X9PPl\n//wbM3tU0qQkgh1Ao4oeJN3EMXZ1CPrw1My2S9op6XshrwsAZRQ9SLqJY+zqECzYzezFkg5J+qi7\n/zrh53vNbMbMZk6dOhXqbQEgVdEZeBPH2NUhSLCbWU/nQv2gu9+R9Bp3v83dp9x9auvWrSHeFgAG\nKjoD37NzUrded5kmJ8ZlkiYnxnXrdZd1qr4uBaixm5lJ+oKkR939s9WHBABh7Nu9Y12NXcqegQ/7\nGLs6hOiK2SXp/ZIeMrPjy9/7hLt/O8C1AdSgC73aIcbYxEHSbWDuPvQ3nZqa8pmZmaG/L4CNnSLS\nuVlsm0oOXRhjE8xs1t2nsl7HlgLAiCnaKRLS4WNz2jV9VJfsv0u7po+mLvxpcowxYEsBYMQ01atd\npKc8ln7ypjBjB0ZMU73aRWbhsfSTN4VgB0ZMUq+2JD37uzO17olSZBYeSz95Uwh2YMSs9Gpv2dxb\n9/35hcVaN7wqMguPpZ+8KXTFACNq1/RRzSXMlicnxnXf/quCvx+dLtXl7Yrh4SkworJKI6F73Ue1\np7wJBDswoi6cGE+csV84MV54V8S8YljV2QXU2IERNegBJX3k3UawAyNq0ANK+si7jVIMMMLSSiOD\nyjRoP2bsADagj7zbmLED2IAOlm4j2AEkaqKDpQvbCXcBwQ6gFepqsRxFBDvQYTHNcAe1WHb1MzWF\nYAc6KrYZLi2W4dAVA3RUbIuI2Ko3HIId6KjYZri0WIZDKQZoUJUaeWyLiGixDIdgBxpStUa+b/eO\nxG1wuzzDZZOwMAh2oIIqM+6qXSDMcJGGYAdKqjrjDlEjZ4aLJEEenprZF83sF2b2cIjrAV1QtSsl\nRBfI4WNz2jV9VJfsv0u7po/WemYpuiNUV8yXJF0d6FpAJ1SdcVftAlm5Y5ibX5Dr+TsGwh1Bgt3d\nvyvp6RDXArqi6oy76oHNsfWxIxxq7EBJIbpSqtTIy9wxxLQFAdINLdjNbK+kvZK0bdu2Yb0tUIuV\ngFxYXNKYmZbcNTmkoFx5b0/5edodQ2xbECDd0ILd3W+TdJskTU1Npf1/EmiF/pntm161Vfc+dkpP\nzS9oYnNPvzq9qLPLr11yV2+TFQ71MrPn/nDuN+iOoa2bbHEXER6lGKBP0sz2P+5/YvXnz5xe3PA7\ni2ddt9x5IncglZ09J4Xziqw7hjZuQcBdRD2CBLuZfVnSH0u6wMyelHSzu38hxLWBYRsUnoPML2wM\n/CLvkWf2nBbCJum+/Vet+17/THhicy/xj1KTWxC09S6i64IEu7u/N8R10F6jdLs8jBls2dlz3v1h\nkmbCvU2m3phpcen5SmjTWxC08S4iBuzuiEyj1i9ddga7ZXMv8ftJi4jKtkrm7X1PmgkvnnW96Pzz\n1rVXXv+6SR04crKxBU5s1VsPgh2ZRq1f+k2v2lr4d3pjppuvuXTD99P+KG7/vXPhulae2XPe3ve0\nGe+vFhZ13/6r9NPpt2vf7h06NDvX6B9stuqtBw9PkWnUbpfvfexU4vfHzHTWfUOXzKDSVNofxf/+\n8dPr2hVN0vWvy9fTnqf3PU/Jpg31bTYyqwfBjkyx7fudJe0P1ll3/XT67UGu1d/v60r/g1JGnsVT\nbfmDzUZm4VGKQaYqt8td3KSqTN037XMW+eMXMlDzlGyob8eLGTsylb1d7mqPctGtAgZ9zqRrmTbO\n2KXwgZo1E47xoA6cQ7AjlzK3y22o4ZZR9A/ZoM+50lvev4r10OxcqUAN2XZKfTteBDtq05Ya7tow\nnNjck/u57pBBQVbkD1nW50y61tQrXtaKOyDq23Ei2FGbNjx07Q/DtSsvQ5WGynzOUboDwvDx8BS1\nqbNHOe9D2aztAUL04w+rFzvtziDpjwpGG8GO2lQ9SCJNkZWweco+VUtDdX3Ofml3ACZ1otsIw0Mp\nBrWqo4ZbpCSRVibpf01VdXzOpK2DD97/RGIPPOUYrMWMHZ1T5KFsUplkrba29yXdlRyanUs9XCPW\nVcAoh2BH5xRZWNNfJtmyuaeJ8V6tJZMQ0u5Kxqx/h5lzWFSEtSjFoHOKLqzpYktf2gx8yV3jvTEW\nFWEgZuzonGE9rGxS2gx85bMW/exd3NoB5Zn78I8fnZqa8pmZmaG/L9C0vCtHk842He+NlfoDFvJa\naJaZzbr7VNbrKMUgaoOCNE/IZr2myBL/IitHqy73XzuuTWZa6pvAsbApbgQ7ojUoSCVlhmxWEBdd\n4l905WjZZwP94+oP9RV00sSLYEe0sk5+ygrZtN//1DdPaM/Oycyg7p/Np/XThw7YvIdx00kTL4Id\n0SqzCdnan6W97pnTizp8bG7g9ZNm88ParjfPHwo6aeJGVwyiNajfPU8v/KDAPXDk5MBrJM2aXSp1\nzmlRaeMaM4u2iwjrEeyI1qDNufJs3DUocJ+aXxh4jUFH4tXdppk2rn/688v10+m36779VxHqkQtS\nijGzqyX9i6QxSf/u7tMhrovuCHkARCh5OksG/WzPzkndcucJzS8sbrj2S8d7A69/4MjJxJr65MT4\n6uEbdeEADVTuYzezMUk/lPRWSU9K+r6k97r7I2m/Qx97XLrYJ12kn3zf1x7Q4tn1/570xkwH3nV5\n7tZGqf3/naD98vaxhyjFvF7Sj9z9J+7+nKSvSLo2wHXREVndJ21TZNvfPTsn9eIXbryxXVzygZ9v\nFFbHor1ClGImJf1szddPSnpDgOuiI9pyBF5eRfvJ509vLMVI2Z+vi3vUIA5De3hqZnvNbMbMZk6d\nOjWst8UQFNltsQ2K/iHq2ucDQgT7nKSL13x90fL31nH329x9yt2ntm7dGuBt0RbDOhoulKJB3bXP\nB4QI9u9LeqWZXWJm50t6j6Q7A1wXHdG1enKZoH7Bec//q7Jlc0+3XneZJLFjIlqpco3d3c+Y2Ycl\nHdG5dscvuvuJyiNDp9RdTw7ZTlmkHTCpu+W3i2c18/jTOjQ7l3ufGGCY2LYXwdTVy56ndbCu9941\nfTSxH30sYcdEaTh96hhdw2x3BAq1EBaV1U5Z53sPOsmoyOuBYSLYEUSdvexZm23d8NUHanvvop0v\nLx3vVX5PoCqCHUGkhe/c/ELlB4tp4Tqxuad9X3ug1tlz0oPWQZ597gwPUdE4gh2VrJylOehJzdz8\ngvb95wOlAy+ti+W3i0sblvqvlTXbznMO6ErHz5j178uYLGtFKjAMBDtKW1vbzrK45PrUN4s3S608\nFF1YXFoN15V2yoXFs6m/l9W+WHRbgbMFmgyos6NpBDtKy3tSz4pnUpbmp+n/w7HkvhrYWR0vWX30\nRZ8JFKm1syIVTSPYUVrdM9Os8N2yOflB5ZbNvczgL7qtQFI5qLfJ1BtbX6JhRSragGBHaUVnphMF\nO0aywvfmay7dEKy9MdPN11yaee2i2wokra498O7LdeBdl3dmxS1GB2eeYoO8i3327d6xYeHQILe8\nIztw10o7AHolfKscKJE09qzZdtrqWoIcbUOwY52kQ5jTlsqvDdZBhzVL0iYrHoB5wrfsVgacMoSY\nEexYp+he5WuD9fCxOX3ijgd1OqFb5S/fsK3wWPbsnNTM40/ry9/7mZbcNWam619XbU+aNh7hB4RG\njR3rVDk0Y8/OST3y6bfpfVduW21NHDPT+67cps/suazwWA4fm9Oh2bnVBUhL7jo0O1e6H77OrQeA\nNmETMKyTtulVqM2tisyY08ayZXNPm88/r/Csu+7PBtSNTcBQSp2HShSdMactfHrm9GKpWXfXjvAD\nyiLYsU6dh2YUXRSUdxl/3g2/OOIOo4KHpx0w7Ad+oQ/NWBl/2gy86Na4Ra6xVpkWxzJ4QIumEewt\nV6T9sI2SDsnolzZjnkzpYy9yjbWG0eLY9f+9EAeCveWKth+2TdZ+MoNmzEkz7N6YSa51uzoWmXXX\nfYRf1//3QhwI9pYbtM/5Jfvvqu1WP1Q5YVCJZDLjumkz7KTvtSU0eUCLNiDYWy5tWb2kdV0hUrhb\n/ZDlhLTx520x7Noy/qxtEIBhoCum5fKc4BPqGLgVIY+5q7N9so1G7fOinZixt1x/OSKtTyTkrX6V\nckJSCefW6y5rbekkNPagQRsQ7B2wthyRtnoy5K1+2XJCWgnn1usuG6mVnXU/oAWyVCrFmNm7zeyE\nmZ01s8xlrqhuGLf6Zd8jZAkHQHlVZ+wPS7pO0r8FGAtyGMatftn3oCMEaIdKwe7uj0qS5Vz6jTCG\ncatf5j3SSjibzHT42Bzb7QJDQlcMCjt8bE67po/qkv13adf00dUNuNI6eJbcK22Py3a7QDGZM3Yz\nu0fSyxN+dJO7fyPvG5nZXkl7JWnbtuKHLiC8MrPgpAekH7v9uGYef3p1z/UbvvrAhn1eqqy+ZDUn\nUExmsLv7W0K8kbvfJuk26dx+7CGuifLKLkJKClmXdPD+JzT1ipdpz85Jfez244m/W7bWTu0eKIZS\nTMeklUGKKtvBkhamLumjtx/Xrumjeul4L/E1ZVsy2W4XKKZqu+M7zexJSW+UdJeZHQkzLCQJWWsu\nOwvOCtO5+QU9+9wZ9Tatf6BepSWT1ZxAMZWC3d2/7u4XufsL3P333X13qIFho5B94mVnwft271BW\nD9TikuvFLzwv2GEddR7+AcSIlacdErLWXPbQiT07JzXz+NM6eP8TqdsbSNL86UUd++SfFB7XoPcl\nyIF8qLF3SNpseqVPvIgqs+DP7LlM//wXV2hywOye+jfQHPMCx4+FMjU15TMzM0N/364bdBrReG+s\nkfJE0piaGgsQOzObdffM7VsoxXTISlCG7hMPMaa8/fCsIAXqR7B3TB194lXlrX9zHigwHNTYO6ir\nfd3s/ggMB8HeQU31dVddHMUKUmA4KMXUqK56chOn9IQoo3AeKDAcBHtNygRhkT8Ew+7rDrERV9ne\neQDFUIqpSdF6ctu3pg1RRmEFKTAczNhrkhWE/bPz08+dafXWtKHKKKwgBerHjL2EPA8RB3WuJM3O\nnzm9mPj6tjxYZCMuoDsI9oLylkwGBWFSmSZNWx4sUkYBuoNSTEF5HyIO6lxJW2DUL++MeFirOSmj\nAN1AsBdU5CFiWhCm1asnxnt60QvOq3xUHas5gdFGsBcU4iFiWtvfLe+4tHAYcx4ogH7U2AsK8RAx\nZL2a1ZwA+jFjL2hQ7byJBUas5gTQj2AvISmUm6p1s5oTQD9KMYE0tXMhbYgA+jFjD2RYte60cg9B\nDmAFM/ZAhrFHetv3kwHQDgR7QWnbCQxjyT0HVQDIo1IpxswOSLpG0nOSfizpr919PsTA2ijPA9I6\nV4DS2gggj6o19rsl3ejuZ8zsHyXdKOnvqw+rnbIWA9Vd66a1EUAelUox7v4ddz+z/OX9ki6qPqT2\nanrGzA6LAPII2RXzQUm3B7xe69Q5Y86zuKmJI/EAdE9msJvZPZJenvCjm9z9G8uvuUnSGUkHB1xn\nr6S9krRt27ZSg21aXYuBiixuorURQBZz92oXMPuApL+V9GZ3P53nd6ampnxmZqbS+zalji1yd00f\nTbwTmJwY1337r6p07bKGtRUwgPzMbNbdp7JeV7Ur5mpJH5f0R3lDvevqmDE3Xbvvx1bAQLdV7WP/\nnKSXSLrbzI6b2ecDjGnkDGNxUxH0ywPdVrUr5g/d/WJ3v2L5nw+FGtgoaVu3S9vuIAAUw8rTFmjb\nRl5tu4MAUAybgLVEW7pdDh+b07O/O7Ph+/TLA91BsGNV/0PTFVs293TzNcWP7QPQDEoxWJX00FSS\nNp9/HqEOdAjBjlU8NAXiQLBjFQ9NgTgQ7FjVtrZLAOXw8BSr2GQMiAPBjnXa0nYJoDxKMQAQGYId\nACJDsANAZAh2AIgMwQ4AkSHYASAyBDsARIZgB4DIEOwAEBmCHQAiM7JbChw+NseeKACiNJLB3n9S\n0Nz8gm684yFJItwBdN5IlmKSTgpaWFzSgSMnGxoRAIQzksHOSUEAYlYp2M3s02b2oJkdN7PvmNmF\noQZWJ04KAhCzqjP2A+7+Wne/QtK3JH0ywJhqx0lBAGJW6eGpu/96zZcvkuTVhjMcnBQEIGaVu2LM\n7B8k/ZWkX0l6U+URDQknBQGIVWYpxszuMbOHE/65VpLc/SZ3v1jSQUkfHnCdvWY2Y2Yzp06dCvcJ\nAADrmHuY6omZbZP0bXd/TdZrp6amfGZmJsj7AsCoMLNZd5/Kel3VrphXrvnyWkmPVbkeAKC6qjX2\naTPbIemspMclfaj6kAAAVVTtirk+1EAAAGGM5MpTAIgZwQ4AkSHYASAyBDsARIZgB4DIEOwAEBmC\nHQAi06mj8TinFACydSbYOacUAPLpTCmGc0oBIJ/OBDvnlAJAPp0Jds4pBYB8OhPsnFMKAPl05uEp\n55QCQD6dCXaJc0oBII/OlGIAAPkQ7AAQGYIdACJDsANAZAh2AIgMwQ4AkTF3H/6bmp2S9PjQ3zic\nCyT9sulB1IzPGAc+YxxWPuMr3H1r1osbCfauM7MZd59qehx14jPGgc8Yh6KfkVIMAESGYAeAyBDs\n5dzW9ACGgM8YBz5jHAp9RmrsABAZZuwAEBmCvQQzO2Bmj5nZg2b2dTObaHpMdTCzd5vZCTM7a2bR\ndB2Y2dVmdtLMfmRm+5seTx3M7Itm9gsze7jpsdTFzC42s3vN7JHl/59+pOkxhWZmLzSz/zWzB5Y/\n46fy/B7BXs7dkl7j7q+V9ENJNzY8nro8LOk6Sd9teiChmNmYpH+V9DZJr5b0XjN7dbOjqsWXJF3d\n9CBqdkbSDe7+aklXSvq7CP+3/J2kq9z9cklXSLrazK7M+iWCvQR3/467n1n+8n5JFzU5nrq4+6Pu\nHttp4a+X9CN3/4m7PyfpK5KubXhMwbn7dyU93fQ46uTuP3f3Hyz/599IelRSVAc2+Dn/t/xlb/mf\nzAejBHt1H5T0X00PArlNSvrZmq+fVGRhMIrMbLuknZK+1+xIwjOzMTM7LukXku5298zP2KkTlIbJ\nzO6R9PKEH93k7t9Yfs1NOnc7eHCYYwspz+cE2szMXizpkKSPuvuvmx5PaO6+JOmK5Wd5Xzez17j7\nwGcnBHsKd3/LoJ+b2Qck/ZmkN3uHe0azPmeE5iRdvObri5a/hw4ys57OhfpBd7+j6fHUyd3nzexe\nnXt2MjDYKcWUYGZXS/q4pHe4++mmx4NCvi/plWZ2iZmdL+k9ku5seEwowcxM0hckPerun216PHUw\ns60rXXdmNi7prZIey/o9gr2cz0l6iaS7zey4mX2+6QHVwczeaWZPSnqjpLvM7EjTY6pq+aH3hyUd\n0bmHbV919xPNjio8M/uypP+RtMPMnjSzv2l6TDXYJen9kq5a/vfwuJn9adODCuwPJN1rZg/q3KTk\nbnf/VtYvsfIUACLDjB0AIkOwA0BkCHYAiAzBDgCRIdgBIDIEOwBEhmAHgMgQ7AAQmf8HGGo7w/4b\nx8sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1182dd990>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The probability of HSIC being its size or larger under the hypothesis that X and Y are independent is 0.0\n"
     ]
    }
   ],
   "source": [
    "X = np.random.randn(100,1)[:, None]\n",
    "Y = X + 0.5*np.random.randn(100,1)[:,None]\n",
    "\n",
    "p_value = HSIC_test(X, Y, gaussian_kernel, 1, gaussian_kernel, 1)\n",
    "\n",
    "plt.scatter(X[:,0], Y[:,0])\n",
    "plt.show()\n",
    "print \"The probability of HSIC being its size or larger under the hypothesis that X and Y are independent is \" + str(p_value)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
